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PREFACE 


This  report  was  prepared  by  Dr.  G.L.  Guymon  and  T.V.  Hromadka  II 
of  the  School  of  Engineering,  University  of  California,  Irvine. 

This  report  covers  work  funded  by  DA  Project  AA7627 19AT42 , Design, 
Construction  and  Operations  Technology  for  Cold  Regions,  Task  A3, 

Facilities  Technology/Cold  Regions,  Work  Unit  006,  Volume  Change  Induced 
by  Freezing  and  Thawing  of  Pavement  Systems.  The  research  consisted  of 
two  components:  1)  the  initial  stages  of  development  of  a two-dimensional 
and  three-dimensional  heat  transport  model  (with  phase  change)  for  freezing 
soils,  and  2)  the  continued  development  of  a one-dimensional  frost  heave 
model  based  upon  solution  of  the  coupled  heat  transport  and  fluid  transport 
problem.  Only  the  first  element  of  research  is  reported  herein.  The 
second  element  of  the  research  is  being  reported  separately  in  the  form  of 
sections  of  another  report. 

It  is  emphasized  that  this  report  covers  the  initial  phases  of  the 
development  of  a multidimensional  heat  transport  model.  While  the 
techniques  used  herein  are  valid,  it  is  expected  that  subsequent  versions 
of  the  model  will  be  prepared  to  increase  efficiency  and  accuracy.  In 
particular,  more  accurate  techniques  of  handling  the  phase  change  problem 
can  be  envisioned  and  these  techniques  will  be  likely  to  reduce  computer 
storage  requirements  and  require  less  computer  time.  Subsequent  improvements 
of  the  model  will  be  reported  upon  in  appropriate  CRREL  reports. 

Cameron  Appel  of  CRREL  technically  reviewed  the  manuscript  of  this 
report. 
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Chapter  1 
INTRODUCTION 


This  report  presents  a model  of  transient  heat  conduction  In  a freezing 
and  thawing  soli.  The  partial  differential  equation  for  transient  heat 
conduction  Is  solved  by  a finite  element  analog  using  a quadratic  weighting 
function  for  the  discretized  spatial  domain.  The  transient  problem  Is  solved 
by  the  Crank-Nlcolson  approximation.  Phase  change  Is  approximated  as  an 
Isothermal  process. 

Both  a two-dimensional  and  three-dimensional  model.  Incorporated  In  the 
same  computer  program,  are  presented.  In  the  latter  case.  It  Is  anticipated 
that  certain  problems  can  only  realistically  be  modeled  as  a three-dimensional 
system.  Examples  of  such  problems  Include:  thaw  degradation  around  roadway 
culverts,  embankment  dams  on  permafrost  where  the  dam  length  Is  short  relative 
to  the  dam  width,  and  thaw  or  freezeback  under  buildings.  However,  In  many 
cases  the  more  economical  two-dimensional  model  may  be  used.  Examples  of  such 
problems  Include:  embankment  dams  of  great  length,  roadway  cross-sections, 
and  long  pipeline  problems. 

This  report  develops  the  basic  equation  of  heat  transport  and  the 
assumptions  and  limitations  upon  which  the  model  Is  based.  The  finite  element 
method  Is  reviewed  and  a complete  derivation  of  the  system  analog  Is  presented. 
Numerous  model  evaluations  were  made  and  are  summarized  herein.  It  Is  empha- 
sized that  the  primary  thrust  of  evaluations  Is  the  numerical  testing  of  the 
model.  Field  verification  has  not  been  attempted  yet.  A user  manual  and 
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computer  listing  of  the  program,  written  In  standard  FORTRAN  IV,  are  pre- 
sented. Although  the  model  was  developed  on  an  IBM  370/155  computer,  there 
should  be  little  difficulty  In  adapting  the  model  to  other  computers  such  as 
the  Dartmouth  System. 

The  main  advantage  of  the  model  presented  here  is  that  it  can  be  readily 
adapted  to  complex  shapes  which  sometimes  is  a problem  with  the  finite 
difference  methods.  The  model  can  accommodate  variable  element  sizes  and 
configurations  using  triangular  shaped  and/or  rectangular  shaped  elements  in 
the  two-dimensional  case  and  tetrahedra  shaped  and/or  brick  shaped  elements 
in  the  three-dimensional  case. 

The  program  has  been  prepared  in  a highly  efficient  manner,  minimizing 
as  much  as  possible  computer  execution  time  and  minimizing  the  storage 

required  for  arrays. 

I 

I 

Future  work  on  the  model  will  require  field  testing  and  verification. 
Additionally,  it  is  desirable  to  couple  a more  sophisticated  boundary  condi- 
tional routine  to  this  model  in  order  to  more  readily  simulate  the  soil  air 
interface.  The  present  model  only  handles  a specified  boundary  condition  or 
I a no  heat  flux  boundary  condition. 
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Chapter  2 
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HEAT  CONDUCTION  EQUATION 

A rigorous  derivation  of  the  heat  transport  equation  can  be  found  in  Bird 
et  al.  (1960).  For  purposes  of  this  report  a more  simplistic,  but  correct, 
derivation  will  be  presented.  Generally,  most  references,  for  example, 

Myers  (1971),  begin  by  making  the  deterministic-continuum  assumption  which 
usually  leads  to  a partial  differential  equation  with  temperature  as  the  state 
variable  and  various  state  parameters,  e.g. , heat  capacity,  thet  arise  out  of 
necessary  mathematical-physical  assumptions. 

The  first  concept  that  needs  to  be  employed  is  that  energy  is  conserved. 
Thus,  by  considering  the  various  rate  processes  Involved  in  a particular 
process  and  by  making  an  energy  balance  on  a control  volume,  the  appropriate 
h6at  equation  is  obtained.  The  various  rate  processes  that  might  be  considered 
are:  conduction,  convection,  radiation,  heat  storage,  and  heat  generation  (e.g. 
latent  heat  effects). 

The  primary  processes  in  a given  soil  system  Include  all  of  these  processes 
if  the  soil  is  freezing  or  thawing.  Moreover,  the  soil  system  Includes  a 
heterogeneous  mixture  of  dissimilar  materials:  mineral  soil,  organic  material, 
air,  water,  and  ice.  Moreover,  the  water  is  often  a dilute  solution  containing 
dissolved  minerals  which  affect  the  system's  thermal  properties.  The  thermo- 
dynamics of  soil  systems  is  treated  by  Edelfsen  and  Anderson  (1943)  among 
others. 

The  derivation  below  will  ignore  radiation  since  this  process  occurs  at 


3 


the  soil  surface.  Radiation  will  be  Included  in  the  system  boundary  conditions. 
Momentarily  we  will  Ignore  energy  generation  in  the  soil  due  to  freezing  or 
thawing.  Thus,  consider  a three-dimensional  elemental  volume  of  material  in 
the  presence  of  a fluid  flux  field.  That  is,  fluid  is  moving  through  the 
elemental  volume. 

The  energy  balance  equation  is 


E 

c 


+ E 


V 


(1) 


where  E is  the  net  rate  of  heat  conduction  into  the  elemental  volume,  E 
c V 

s the  net  rate  of  heat  convection  into  the  elemental  volume,  and  E^  is  the  total 
ae  heat  energy  is  stored  in  the  elemental  volume.  Consider  the  x-dlrectlon 
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= net  heat  conduction  and  convection  for 
the  x-dlrection 


where  AxAyAz  is  the  volume  of  the  element,  T^  Is  a reference  temperature, 
c is  the  volumetric  heat  capacity  of  the  fluid  and  v is  the  fluid  flux 

w X 

in  the  x-directlon.  The  variabJ'*  q may  be  replaced  with  Fourier's  Law;  l.e.. 
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where  k.  Is  the  thermal  conductivity  of  the  entire  mass  of  material  in  the 


elemental  volume.  Substituting  this  law  into  the  above  yields 


9x 


(k  1^)  - (c  V T) 

\ X dx  / dx  w X 


for  the  rate  of  conduction  and  convection  in  the  x-dlrection.  Now  considering 
the  y-  and  x-directlons  in  turn,  similar  expressions  are  added  together  and 
equated  to  the  rate  of  heat  energy  accumulations  in  the  elemental  volume.  Thus, 

I;  {k  f ) S)  - i M 
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where  c^  is  the  volumetric  heat  capacity  of  the  mixture  and  0 is  time. 

Now  looking  at  the  convection  term  and  expanding  by  the  chain  rule 

assuming  an  incompressible  fluid  (i.e.  c = constant) 

w 


/ 9 (v  T)  9 (v  T)  9 (v  T)  \ 
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c (v  -p+v  |£+v  + + ^ + 
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However,  the  last  term  is  equal  to  zero  since  it  is  none  ether  than  the 
continuity  expression  for  an  incompressible  fluid.  Thus,  the  final  result  is 


/ at  ^ 3T  3T  \ 3T 

- % (''x  9^  ^ ^ ''z  ^ j = "a  ^ 


A major  assumption  that  will  be  incorporated  in  the  model  developed  herein 
is  that  the  model  will  be  developed  for  a system  in  which  fluid  flow  is 
negligible.  Therefore,  the  convective  terms  will  be  Ignored  and  the  following 
equation  will  be  solved  in  general: 


/k 

1 + ^1 

fk 

(v 

9T 

3x 

( X 9x  j 

9y  1 

y 3y/  9z 

( z 9z/ 

- *^3  90 

(3) 

This  equation  will  be  solved  in  both  three-  and  two-dimensions.  The  two- 
dimensional  form  is  derived  from  Equation  3 by  assuming  9T/3Z  = 0 and  that 

in  the  z-directlon  the  system  consists  of  a slab  of  uniform  thickness.  Thus, 


i_  /k  3I\  + i_  /v  i!I\  = il 

9x  \ X 3x  / 3y  \ y Sy  / '^a  99 


Both  equations  are  of  the  parabolic  type.  Boundary  conditions  to  be  considered 
are  specified  conditions  and  no  heat  conduction  conditions.  Initial  conditions 
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are  also  required.  Latent  heat  generation  will  be  approximated  as  an  Iso- 
thermal process  as  will  be  described  subsequently. 


Chapter  3 


FINITE  F.LEMENT  METHOD 


Introduction 

The  finite-element  method  has  been  routinely  used  in  the  structural, 
mechanical,  and  aerospace  engineering  fields  since  the  use  of  modern  digital 
computers  became  widespread.  The  finite-element  literature  relative  to  these 
fields  is  extensive.  Following  the  mid-1960's,  investigators  concerned  with 
general  field  problems  became  Interested  in  this  powerful  numerical  tool, 
and  the  volume  of  published  articles  dealing  with  the  numerical  solution 
of  general  problems  by  the  finite-element  method  has  Increased  markedly. 

Several  recent  texts  present  an  excellent  treatment  of  the  finite-element 
method:  Meyers  (1971),  Desal  and  Abel  (1972),  Zienklewlcz  et  al.  (1971), 

Huebner  (1975),  and  Sergerllnd  (1976).  The  latter  two  texts  are  especially 
recommended  for  the  beginner. 

The  finite-element  method  is  Ideally  suited  to  deal  with  complex 
geometries,  anisotropy,  and  heterogeneity  which  are  characteristic  of  most 
practical  problems.  For  certain  classes  of  problems,  the  finite-element 
method  may  be  a more  efficient  numerical  technique  than  the  traditional  finite- 
difference  methods.  That  is,  the  finite-element  computer  program  may  require 
less  execution  time  than  the  finite-difference  computer  program  for  a speci- 
fied level  of  precision.  In  particular,  time-dependent  problems  are  efficiently 
solved  by  the  finite-element  method.  The  finite-difference  methods  often  pre- 
sent stability  and  truncation  error  difficulties.  Finite-element  methods  can 
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readily  use  arbitrary  mesh  spacing  and  can  easily  handle  complicated  boundary 
conditions  that  usually  require  lengthy  programming  by  finite-difference 
methods.  The  finite-element  program,  once  written  for  a class  of  problems, 
is  very  effective  since  it  can  be  used  to  solve  similar  problems  for  any  geo- 
metry and  any  mesh  configuration  desired.  In  other  words,  the  finite-element 
program  is  completely  general  for  the  class  of  problems  it  was  designed  for. 
Finite-difference  programs  are  often  special  purpose  since  they  apply  only  to 

a specific  geometry  with  a specific  mesh  spacing.  For  simple  geometries 
and  steady-state  conditions,  the  finite-element  methods  appear  to  offer  little 
advantage  over  the  more  familiar  relaxation  and  iteration  methods  commonly 
applied  to  Laplace-like  problems.  Unfortunately,  research  on  the  various  nu- 
merical methods  has  not  advanced  to  the  point  where  definitive  criteria  can 
be  stated  for  the  selection  of  the  best  method  to  use  in  a particular  case. 
Indeed,  there  is  still  a good  deal  of  controversy  over  the  comparative  advan- 
tages and  disadvantages  of  the  various  methods. 

There  are  two  general  approaches  to  the  finite-element  method:  (i)  the 
direct  approach  which  Involves  writing  a set  of  system  matrices  by  visualizing 
the  physical  linkages  of  a system  and  (il)  the  variational  approach  which 
Involves  developing  a variational  principle  or  applying  the  Galerkln  tech- 
nique. The  direct  approach  is  primarily  applicable  to  structural  engineering. 
The  second  approach  consists  of  formulating  a variational  principle  or  for- 
mulating a governing  partial  differential  equation  which  can  be  converted 
to  an  equivalent  variational  problem.  It  is  the  latter  method,  a mathematical 
abstraction,  that  is  applicable  to  general  field  problems  such  as  heat 
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transport. 


For  the  problem  considered  In  this  report,  the  Galerkin  method  and  the 
variational  principle  method  lead  to  identical  results.  The  variational 
principle  method  will  be  used  to  solve  the  problem.  An  extremum  problem  re- 
places the  given  partial  differential  equation,  and  a functional  is  found 
such  that  the  extremum  function  also  satisfies  the  given  differential 
equation  and  its  auxiliary’  conditions.  That  is,  given  the  following 
functional  in  two-dimensional  Cartesian  coordinates 

//f  ^x,y,T,  1^,  l^jdxdy;  in  R 

find  a function  T in  R such  that  the  functional  is  a minimum.  If  such  a 
function  exists,  application  of  the  fundamental  lenma  of  Integral  calculus 
yields  the  so-called  Euler  equation;  l.e., 

|^^3T/8x)  + l^l^aT/By)  "11=0 

It  is  this  equation  that  is  given  to  start  with  for  field  problems.  In  order 
to  develop  the  finite-element  method,  we  resort  to  a mathematical  abstraction 
We  convert  the  governing  partial-differential  equation  describing  our  prob- 
lem to  an  entirely  new  problem,  a variational  problem. 

In  general,  the  finite-element  method  will  be  applied  to  the  linear  two- 
dimensional  parabolic  partial-differential  equation  for  nonsteady  heat 
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transport  In  an  incompressible  porous  medium.  For  generality,  a source 
term  will  be  included  and  the  medium  will  be  considered  anisotropic  and 
heterogeneous.  The  method  can  readily  be  specialized  to  the  one-dlmenslonal 
case  or  steady-state  case  and  can  readily  be  extended  to  the  more  general  three- 
dimensional  case.  The  Cartesian  coordinate  system  will  be  used,  assuming  the 
reader  can  extend  the  derivations  here  to  other  orthogonal  coordinate  systems. 

In  essence,  a specific  problem  is  chosen  as  an  example  around  which  to  develop 
the  salient  aspects  of  the  finite-element  technique. 

A Variational  Principle 

Many  problems  of  heat  flow  in  porous  materials  are  represented  by  the  par- 
tial-differential equation 


where  T is  a continuous,  single-valued  function  (the  unknown  state  variable); 

k and  k are  the  thermal  conductivities  In  the  x and  y directions,  respectively; 
X y 

c is  the  heat  capacity;  Q Is  a generalized  source  term;  and  0 Is  time.  An 
equivalent  variational  functional  is 


ff\>m  ^ - ’) 


dxdy 


11 


where  9T/30  is  assumed  invariant  or  replaced  by  a finite  difference  analog, 
such  as  the  Crank-Nicolson  scheme.  However,  9T/36  will  be  considered  as 
invariant  and  the  time-domain  problem  will  be  dealt  with  later. 

Before  continuing  with  the  variational  procedure,  the  problem  must  be  more 
carefully  specified.  This  consists  of  considering  typical  boundary  conditions, 
and  interface  conditions.  Interface  conditions  are  required  by  the  variational 
technique  that  is  to  be  used  for  mathematical  convenience.  It  is  much  too 
restrictive  to  apply  the  variational  functional  to  the  entire  domain  of  interest, 
since  it  would  be  required  that  the  first-order  space  partlals  exist  through- 
out the  domain.  It  is  more  convenient,  and  incidentally  more  practical  for 
field  problems,  to  consider  it  as  applying  to  particular  subregions  of  the 
domain,  where  subregions  are  separated  by  an  interface,  for  example,  abrupt 
changes  in  material  properties. 

Consider  the  connected  domain  R shown  in  Fig.  1,  with  boundaries  and 
interfaces  as  shown.  In  each  subregion,  R^,  the  partial-differential  equation  ; 

holds  and  on  exterior  boundary  surfaces  the  following  boundary  conditions  hold:  | 

9T/3n  = 0 on  r ,t  > 0 

n I 

i 
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on  are 


T ♦ - T - 

'n  'n 


t > 0 


n dn  4 
' n 


k ^ 
% 3n 


t > 0 


where  n Is  the  normal  direction  to  the  surface  and  s is  the  direction  along 
the  surface.  Additionally,  initial  conditions  are  required  throughout  R l.e. 


The  above  auxiliary  conditions  are  applicable  to  a wide  variety  of  situations 
encountered  in  the  field.  These  conditions  are  chosen  as  an  example  and  are 
not  the  only  auxiliary  conditions  that  can  be  dealt  with  by  the  finite-element 
method. 

It  is  more  convenient  to  rewrite  the  variational  prlniclple  in  the 
following  equivalent  form; 


m = 1 


where 


F 

m 


/ 3T  3T\ 

^x.y.T.g^.gy  I 


There  are  M subregions  In  R,  and  m refers  to  a particular  subregion.  Within 
each  subregion  consider  the  parameters  k^,  k^,  pc  and  Q constant.  However, 
these  parameters  may  vary  from  subregion  to  subregion.  As  before,  9T/96  Is 
considered  Invariant.  It  Is  relatively  easy  to  show  by  taking  a small  variation 
of  y l.e.  9x  for  all  admissible  states  of  the  variable  T that  the  above 
approximates  the  governing  partial-differential  equation,  the  boundary 
conditions,  and  the  Interface  conditions.  Admissible  states  of  T are  defined 
as  (1)  T is  continuous  throughout  R,  (11)  the  first  derivatives  of  T are  con- 
tinuous In  R^,  and  (111)  T is  equal  to  the  specified  boundary  conditions  on 

r . 

s 

Finite  Element  Representation  of  Region 

The  above  variational  principle  is  formally  solved  by  the  Ritz  method. 

The  Ritz  method  consists  of  selecting  a trial  sequence  of  functions  that  are 
substituted  Into  the  last  equation.  Such  a sequence  of  functions  Is  obtained 
by  dividing  R Into  an  arbitrary  number  of  finite  elements  that  completely 
cover  the  domain  R as  Illustrated  In  the  two-dimension  connected  domain  shown 
In  Fig.  2.  Notice  that  the  curved  boundaries  are  approximately  modeled  by 
straight-line  segments.  The  triangular  element  is  used  here  although  quadri- 
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t ! (x,  4-Xj  4X,,)/  3 

y-  (Vi^y,  ♦yk)/3 

Fig.  3.  General  triangular  element  with  global  coordinates  (x,y)  and 
t local  coordinates  (C*n) 
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lateral  elements  are  more  efficient  from  the  standpoint  of  computer  storage 
allocations  and  run  times.  Within  each  element  the  unknown  potential  states 
are  approximated  by  a polynomial;  for  example,  a linear  one  Is  used  here; 

T * Bi  + B2X  + B 3y 

Admissibility  requirements  are  met  If,  for  a general  element  as  shown  In 
Fig.  3,  the  polynomial  Is  forced  to  pass  through  the  same  value  at  node  points 
which  are  locally  designated  by  1,  j,  and  k.  This  Is  simply  accomplished  by 
writing  the  following  three  equations  for  the  potential  function  at  each  node 
point: 

= Bi  + B2X^  + Biy^ 

Tj  = Bj  + B2X^  + Bsyj 

\ 

and  applying  Cramer's  rule  to  solve  for  Bi,  B2 , B3.  The  result  for  the  m-th 
element  Is  written  In  the  following  compact  matrix  notation 

T™  = [ A ] {T}™ 
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p 


I 


The  braces  indicate  the  following  column  vector  for  the  potential  at  each 


The  brackets  Indicate  a row  matrix  for  the  m-th  element  of  the  form 


[A  1 = ( A^,  A ] 


The  form  of  A coefficients  are  given  in  detail  in  Zienklewlcz  (1972)  and 
n 


algebraic  relations  are  given  for  the  calculation  of  matrix  element  coeffi- 


cients. These  coefficients  are  functions  of  x and  y,  the  coordinate  position 


of  the  triangle  nodes,  and  the  area  of  the  triangular  element.  In  general. 


braces  will  indicate  a column  matrix  and  brackets  will  indicate  a square 


matrix.  Values  of  the  potential  state,  T,  are  now  defined  in  a unique  and 


continuous  manner  throughout  R. 


To  illustrate  the  finite-element  technique,  a simple  polynomial  is  used. 


Greater  precision  can  be  obtained  for  a particular  element  by  using  a higher 


order  polynomial. 


Minimization  of  Variational  Functional 


The  variational  functional  is  minimized  with  respect  to  the  potential 
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state  T (n  » 1,  ....  N)  at  each  node  point  (where  there  are  N nodes)  by 
n 


evaluating  each  differential  3x/3T^  for  eacTi  element  and  equating  all  such 


contributions  to  zero.  For  example,  the  contribution  over  the  m-th  element 
for  the  1-th  node  Is 


lx 

3T 


3F 


1 ■ 


dxdy 


where 


3F 


3T 


m , 3T  3 /3T\^,  3T  3 /3T\^/„3T  -\3T 

^ “ ‘‘x  3x  3T^  \3x)  S 9y  \3y)  \ 30  " 3T^ 


The  polynomial  shape  function  is  differentiated  with  respect  to  time  to  get 


3T/30  = [A]  {31/30} 


This  equation  together  with  the  polynomial  shape  function  is  substituted  into 
the  above  to  give 


3F 


9t7  ■ “x  H ‘ > w”)  * N 


('■I''!  mr 


Since  [A  1 Is  a function  of  x and  y,  constants  x , y and  the  constant  element 

n n 

area  A , differentiations  with  respect  to  x and  y are  readily  performed, 
tn  ' 

Also  differentiation  with  respect  to  is  easily  performed  since  only 

{T}“  is  a function  of  the  node  state  variables.  Assuming  parameters  are 

constant,  integration  yields  the  desired  result. 

A similar  operation  is  carried  out  for  9x/9T  and  9x/9T  , and  the 

J K 

results  are  combined  in  the  following  compact  matrix  notation 

{9x/9T}”  = Is]  {T}“  + [p]  {9T/90}"'  - {r} 


The  [s  I and  [p  j matrices,  defined  in  numerous  other  references,  are  square, 
symmetrical  3x3  matrices  that  are  functions  only  of  the  global  coordinate 
position  of  the  nodes  on  the  m-th  triangular  element,  the  parameters  for  the 
m-th  triangular  element,  and  the  area  of  the  m-th  triangular  element.  The  {r} 
matrix  is  a function  of  the  source  Q and  the  area  A™. 

The  above  equation  is  general  in  that  it  applies  to  any  interior  element 
of  R provided  there  are  no  specified  boundary  nodes  on  the  element.  For  such 
cases  where  there  are  specified  boundary  conditions,  the  differential  of  x 
with  respect  to  that  node  is  meaningless  since  there  is  no  variation.  This 
difficulty  will  be  resolved  later. 

All  element  system  matrices  are  combined  in  accordance  with  the  following 
equation 

'4^,  ^ 

tn  = i 


which  yields  the  following  system  matrix  expression 


I S*1  {t}  + 1 p*l  {3T/30}-<R*> 


The  brackets  indicate  a square  symmetric  matrix  N X N and  the  braces  indicate 
and  N column  matrix  of  the  unknown  potential  states  including  the  specified 
boundary  nodes.  Essentially,  the  above  is  a formal  statement  that  all  element 
contributions  to  a particular  node  are  added  together  to  form  the  equation 
for  that  node.  The  equation  for  a particular  node  then  appears  in  the  matrix 
in  an  ordered  manner. 


The  above  cannot  yet  be  set  to  zero  since  the  prescribed  boundary  conditions 
have  not  been  considered.  Although  the  boundary  conditions  could  have  been 
considered  at  the  element  level,  they  will  be  handled  at  the  system  level  here. 

A scheme  for  handling  the  specified  boundary  conditions  at  the  system  level  is 
given  in  the  subsequent  chapters.  The  starred  matrices  of  the  above  are 
reformed  by  eliminating  equations  associated  with  boundary  condition  nodes  to 
yield  matrices  [S  1 , (P  I , and  {R},  and  the  previous  equation  Is  equated  to 
zero 


2 


source  term  and  the  given  geometric  boundary  conditions.  The  natural  boundary 
conditions  need  not  be  considered  since  they  are  automatically  taken  care  of  i 

by  the  variational  principle. 

The  finite-element  solution  is  formally  complete  once  the  above  has  been 
developed.  This  equation  is  a system  of  linear  ordinary  differential  equations 
which  can  be  readily  solved  by  a variety  of  standard  methods.  That  is,  the 
potential  state  at  each  unknown  node  point  is  solved  for,  and  it  is  assumed 
that  potential  states  vary  in  a linear  manner  between  nodes.  For  small 
problems  where  the  dimensions  on  the  matrices  are  relatively  small,  formal 
integration  can  be  used  to  arrive  at  the  solution.  However,  for  larger 
problems,  numerical  differentiation  techniques  are  preferable. 


i'  “ " ^ /v^r 


Chapter  4 

DERIVATION  OF  THREE-DIMENSIONAL  TRANSIENT 
HEAT  CONDUCTION  MODEL  WITH  QUADRATIC 
SHAPE  FUNCTION 


Discussion 

This  chapter  examines  a three-dimensional  transient  heat  problem.  The  only 
boundary  conditions  that  will  be  considered  are  (1)  given  surface  boundary 
temperatures  and  (2)  zero  heat  flux  on  the  surface  boundary.  The  solution  will 
be  approximated  by  the  finite-element  method,  using  a variational  statement 
equivalent  to  the  governing  partial-differential  equation.  The  subject  volume 
will  be  discretized  into  three-dimensional  elements,  and  a quadratic  interpolating 
shape  function  assumed  for  the  field  variable  (temperature)  in  each  element.  The 
resulting  information  of  the  problem  is  a finite  set  of  linear  simultaneous 
equations,  the  variables  of  which  are  the  values  of  the  field  variable  at  specific 
interior  nodal  points.  Time  advancement  of  the  solution  vector  at  nodal  points 
is  by  use  of  the  Crank-Nlcolson  method. 

The  computer  program  developed  in  this  report  is  constructed  to  provide  easy 
access  to  the  program  schemes.  Modifications  can  be  made,  and  further  sophis- 
tications incorporated  without  a major  overhaul  of  the  entire  program.  Program 
schemes  are  separated  by  comment  statements  explaining  or  defining  the  compu- 
tation algorlthlms  and  program  variables,  as  they  are  derived  in  this  report. 
Hence,  individual  processes  are  explained  in  detail  by  referring  to  this  report 


while  examining  the  program  printout. 


Fo rmulation 


The  three-dimensional  volume  will  be  discretized  into  tetrahedra  and/or 
"brick"  elements  as  shovm  In  Fig.  4.  The  elements  and  nodal  points  are 
numbered  to  reduce  ultimate  matrix  bandwidth.  The  element  data  consisting 
of  nodal  numbers  and  element  parameters  are  read  into  the  computer  element 
by  element. 

The  model  will  further  subdivide  the  three-dimensional  brick  elements  Into 
five  tetrahedra  elements.  Due  to  the  quadratic  interpolating  model,  there  are 
10  nodes  per  tetrahedron  and  26  nodes  per  brick,  hence  the  program  must  coordinate 
the  bricks'  nodal  numbering  to  the  five  subsequent  tetrahedra  nodal  numbering 
schemes.  The  set  of  five  tetrahedra  resulting  from  the  subdivision  of  a brick 
element  Is  composed  of  four  "corner"  tetrahedra  of  equal  volumes,  each  being 
one-sixth  of  the  brick's  volume,  and  an  interior  tetrahedron  having  a volume 
of  one-third  of  the  brick's  volume.  Each  of  the  five  tetrahedra  will  assume 
the  same  intrinsic  properties  as  ascribed  to  the  original  brick  element. 

The  model  processes  the  tetrahedron  elements,  with  parameters,  and  con- 
structs tetrahedron  element  conduction  and  capacitance  matrices.  These  element 
matrices,  each  being  a symmetric  10  x 10  unbanded  matrix,  are  then  Incorporated 
into  the  overall  global  conduction  and  global  capacitance  matrices,  each  global 
matrix  being  formulated  in  "banded'-'  form  utilizing  the  Inherent  symmetry  of 
the  systems.  The  element  matrices  are  directly  calculated  term  by  term  using 
the  derivations  Included  In  this  report. 

Finally,  once  all  the  brick's  element  matrices  are  combined  Into  the  global 
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Sub  1:  1,20,18.24,11,10,17,19,25,26  1,6, 5, 8 

Sub  2:  5.20,24,22,13,15,1^,26,23,21,  3.6,8,? 

Sub  3:  5.3,1,20,4,9,13,2,11,12  3. 2, 1,6 

Sub  4:  5.1.7.24,9.6,15.8,16.17  3. 1.4. 8 

Sub  5:  5,1,24,20,9.15.13.17,26.11  3. 1.8, 6 

Fig.  4.  Brick  Division  into  Five  Tetrahedron 
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capacitance  and  global  conduction  matrix  system,  boundary  values  are  Inserted, 
the  Initial  condition  Is  Inserted,  and  the  time  progression  scheme  Is  Initiated, 
using  the  Crank-Nlcolson  method  to  move  the  solution  vector  forward  with  time 
Increments  of  A0  . 

Method 

The  volume  will  be  assumed  to  be  In  the  first  quadrant  of  the  three-dimensional 
coordinate  system.  Each  node  Is  numbered  in  a scheme  to  provide  a minimum 
bandwidth  In  the  banded  global  matrices.  Element  data  Is  prepared  by  two 
methods . 

Tetrahedron  element  - the  node  numbering  Is  read  Into  the  computer  In  the 
sequence  shown  in  Fig.  5.  Any  tetrahedron  corner  may  be  used  as  the  first  node 
entry;  however,  the  remaining  nodes  must  follow  the  Illustrated  sequence  to 
ensure  compatibility  of  the  local  coordinate  system  to  the  midside  node  numbering. 

i 

I Each  comer's  global  coordinates  must  be  read  in  as  well  as  the  element's 

parameters . 

! Brick  element  - the  nodal  sequence  of  Input  must  follow  the  pattern  shown 

In  Fig.  6.  It  is  Important  that  the  first  nodal  entry  be  as  shown  In  the  figure. 
Only  the  coordinates  of  the  first  entry  are  read  In;  the  other  corner  coordinates 
are  calculated  In  the  program  using  the  brick  dimensional  data.  From  the  brick's 
nodal  sequence,  five  sequences  are  formulated  In  the  program  for  the  five  ultimate 
tetrahedron  elements;  hence  the  nodal  numbering  Input  must  follow  the  shown  pattern. 

The  tetradedra  utilize  a local  volume  coordinate  system  as  discussed  in  Desai 
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and  Abel.  Utilizing  Fig.  5,  the  Cartesian  coordinates  are  related  to  the 
local  coordinate  system  by  the  relation 


where  each  v Is  defined  to  be  the  volume  of  that  tetrahedron  with  vertices 
at  point  (x,y,z)  and  the  three  nodes  other  than  node  "1"  , 1 = 1,  2,  3,  4 ; 

and  V equals  the  volume  of  the  overall  tetrahedron. 

Define  a^^  to  be  the  cofactor  of  In  the  above  determinant.  Also 

define  b^  and  c^  to  be  the  cofactors  of  y^  and  z^  respectively.  Then 


where  integration  over  the  entire  volume  v simplifies  to 


I 


I 

r 

I 

i 

! 


/ 


q r s 

4 h H 


dv 


V 


p ! q I r ! s ! 6v 
(p  + q + r + s+3)l 


The  cofactors  mentioned  above  are  derived  as  follows: 


1 

1 

1 

ai  = - 

yz 

y3 

y- 

; b2  = - 

Z2 

Z3 

Z4 

hence , 


1 

1 

1 

1 

1 

X3 

X4 

C3  = - 

Xl 

X2 

X4 

Z3 

Z4 

yi 

y2 

y4 

;etc. 


= “ (yszi.  - zayi.  - y2Z^  + Z2y4  + y2Z3  - Z2y3) 
32  = + (y3Z4  - Z3y4  - yiZ4  + ziy4  + yiZ3  - ziya) 


C4  = t (X2y3  - y2X3  - xiya  + yixs  + xiy2  - yiX2i 

Note  that  once  the  cofactors  have  been  determined,  the  volume  of  the  tetrahedron 
can  be  calculated  by  the  expression 

V * xjai  + X2a2  + xaas  + X4a4 
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However,  when  using  a "brick"  element,  the  volumes  of  the  five  resulting 
tetrahedra  are,  essentially,  four  corner  tetrahedra  at  one-sixth  the  brick's 
volume,  and  one  interior  tetrahedron  at  one- third  the  brick's  volume.  Thus, 
using  a "brick"  element  system  easily  solves  the  volume  calculations,  as 
compared  to  solving  for  volumes  by  evaluation  determinants  for  each  tetrahedron 
element. 

Boundary  values  are  Incorporated  as  discussed  in  Myers  (1971) , 

The  global  conduction  and  global  capacitance  matrix  system  are  constructed 
as  the  tetrahedron-element  contributions  are  determined,  A "brick"  element 
matrix  system  formulation  is  not  needed. 


Derivation  of  element  "XK"  and  "XC"  matrices 

The  governing  partial-differential  equation  in  a three-dimensional  transient 
conduction  heat  problem  on  a volume  v , with  boundary  surface  v^  is 


k 

x 


lit 

9x^ 


+ 


pc 


it 

90 


with  boundary  conditions  t (0)  = t^  on  the  surface  S , a subset  of 


and 


dt 

9n 


= 0 , The  initial  condition  is  t (0  = 0)  = t 


1 • 


V 

s 


In  this  relation,  the  units  of  measurement  are: 

k , k , k = thermal  conductivities  in  z,  y,  x-dlrectlon 
z y X 

Btu/hr  - ft  - °F 
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t = temperature,  °F 
0 = time,  hours 

3 

p = density,  Ibm/ft 
c = specific  heat,  Btu/lbm  - °F 

A variational  statement  for  this  three-dimensional,  transient  conduction  problem 
is  the  minimization  of  a volumetric  integral  over  the  volume  v as  shown  below 

I - [k,  (U)  <« 

V 

Equation  6 must  be  minimized  for  every  Instant  in  time  while  satisfying  the 
boundary  and  initial  conditions  stated  in  Equation  1.  The  volumetric  integral 
expressed  in  Equation  6 can  be  equated  into  the  sum  of  four  Integrals 

■ ■ \ {ry ) 

V V 

2 

“vH/2/pc||ldv  (7) 

V V 

or  In  a different  notation 


V 


The  first  step  Is  to  divide  the  volume  v into  "m"  tetrahedra.  Then  the 
integral  in  Equation  7 is  equal  to  the  sums  of  the  integrals  over  each  element 


m m 


e =»  1 e = 1 


■Z)  C "x  (t)  ■*’'  “y  (t) 

e = 1 v*  V* 

+ \ (f)  Ir 

V*-  V*- 

To  evaluate  the  Integrals  stated  above,  the  field  variable  "t"  is  approximated 
by  a quadratic  polynomial.  This  approximation  will  assume  "exact"  values  of 
the  field  variable  at  specified  points  within  the  element,  these  points  are 
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the  element's  nodal  points  as  discussed  In  Desal  and  Abel  (1971)  and  Meyers  (1971), 
A typical  tetrahedron  element  v^“  with  corner  coordinates  P^,  P^,  P^  and 
P^  Is  related  to  the  local  volume  coordinates  defined  as 


^1  ^(e)  • 


i = 1.  2,  3,  4, 


In  the  quadratic  approximation  of  the  field  variable 


t = [Ni,  N2 , ...  , Nio  Ht^} 


where  the  t^  are  the  values  of  temperature  at  the  respective  node  numbers,  and 

[n]  = 

L'  ‘ V \ ‘ \ ’ 3/  ' ■*  ‘^  / \ 1 2/  \ • 

/y,.v/.  V/  v/  \”1 

(10) 


/ ^ \ 

/ . ^ \ 

/ 2 \ 

/ 2 , 

/ V 

/ 

PS  - s) 

^ - \) 

p.-s) 

K-^) 

r.s) 

f4L  L 

\ 1 3 

(‘^^)  Ks)  (‘s\)  (w) 
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Therefore,  for  example, 


3t 

3x 


jv)£  ‘i 


3t 


9L. 


* ““1 

i = 1 


ai 


6v 


(e) 


a? 


6v 


(e) 
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6v 
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6v 


(e) 


(4Li  - 1)  ti  + 4L2t‘  + 4L3t6  + 4L4t'; 


(4L2  - 1)  t2  + 4Lit5  + 4L3t0  + 4L4tio 


(4L3  - 1)  t3  + 4Litf,  + 4L2t8  + 4L4t9 


(4L4  -1)  t4  + 4Lit7  + 4L3t9  + 4L2tio 


. 9t  9t  1 _ j ^ 9t  T j t-  u 

where  , -s—  are  related  to  except  the  a are  replaced  by  b.  , c. 

dy  dz  dx  1 11 

respectively;  (and  is  replaced  by  , k^  respectively). 

Substituting  Equation  10  Into  Equation  9 transfoims  the  transient  heat 
conduction  relation  within  an  element  "e"  Into  a function  of  nodal  point 
temperature  values.  This  process  Is  repeated  for  each  element  "e"  of  the  dis- 
cretized volume  V.  The  resulting  relations  are  combined  and  minimized  with 
respect  to  each  nodal  variable,  producing  a system  of  linear  equations.  The 
boundary  conditions  and  Initial  conditions  are  Inserted.  Then  values  for  nodal 
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temperatures  are  computed  at  specified  time  step  intervals  by  the  Crank- 
Nicolson  time  advancement  routine.  Mathematical  details  are  contained  in 
Appendix  A. 


i 
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Chapter  5 


DERIVATION  OF  TWO-DIMENS lONAI.  TRANSIENT  HEAT  CONDUCTION 
MODEL  WITH  QUADRATIC  SHAPE  FUNCTION 

Discussion 

This  chapter  examines  the  two-dimension  transient  heat  conduction  problem. 

As  mentioned  In  Chapter  4,  only  zero  heat  flux  on  the  surface  boundary  and 
specified  nodal  boundary  temperatures  will  be  considered  as  boundary  conditions. 

Again,  the  variational  principle  as  applied  to  the  three-dimensional  case  will 

i 

be  utilized  In  two-dimensional  elements,  and  a quadratic  shape  function 
assumed  for  the  field  variable  (temperature)  In  each  element.  The  resulting 

formulation  of  the  problem  is  a finite  set  of  linear  simultaneous  equations,  I 

1 

the  variables  of  which  are  the  values  of  the  field  variable  at  specified  nodal  ! 

points.  Time  advancement  of  the  solution  vector  at  said  nodal  points  is  by  j 

1 

use  of  the  Crank-Nicolson  method.  j 

Formulation 

The  two-dimensional  volume  will  be  discretized  into  triangles  and/or 
rectangular  elements.  The  elements  and  the  nodal  points  are  numbered  to  reduce 
ultimate  matrix  bandwidth.  The  element  data  is  read  into  the  computer  element 

by  element  as  to  nodal  numbering  and  thermal  parameters. 

I 

[ The  program  subdivides  the  two-dimensional  rectangular  elements  into  two  triangles 

I of  equal  volume.  The  quadratic  shape  function  employs  nine  nodes  per  rectangle, 

and  six  nodes  per  triangle,  hence  the  program  must  coordinate  the  rectangle's  nodal 

i 

I 

t 

i 
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numbering  to  the  two  subsequent  triangle  numbering  schemes.  The  two  resulting 


triangles  will  assume  the  same  intrinsic  properties  as  ascribed  to  the  two- 
dimensional  rectangular  element. 

Method 

The  area  will  be  assumed  to  be  in  the  first  quadrant  of  the  Cartesian 
coordinate  system.  Each  node  Is  numbered  in  a scheme  to  provide  a minimum 
bandwidth  in  the  ultimate  global  matrix  system.  Element  data  is  prepared 
by  two  methods: 

Triangle  element  - the  node  numbering  is  read  into  the  computer  in  the 
sequence  shown  in  Fig.  7.  Any  triangle  corner  may  be  used  as  the  first 
node  entry,  however,  the  remaining  nodes  must  follow  the  Illustrated 
sequence  to  ensure  compatibility  of  the  local  coordinate  system  to  the 
midside  node  numbering.  Each  corner's  global  coordinates  must  be  read 
in;  as  well  as  the  element's  parameters. 

Rectangle  element  - the  nodal  sequence  of  input  must  follow  the  pattern  shown 
in  Fig.  8.  It  is  Important  that  the  first  nodal  entry  be  as  shown  in 
Fig.  8.  Only  the  coordinates  of  said  first  entry  are  read  in. 

The  triangle  will  utilize  the  local  area  coordinate  system  as  discussed 
in  Desal  and  Abel  (1972).  Utilizing  Fig.  7,  the  Cartesian  coordinates  are 
related  to  the  local  coordinate  system  by  the  relation: 


37 


1 1 l"|  /LA 

Xl  X2  X3  \L2/  , Lj^ 

y 1 72  7 3 \ L3  / 


L = — ; 1 = 1,  2,  3 

A 


where  each  Is  defined  to  be  the  area  of  the  triangles  shown  In  Fig.  9, 


and  A 

= 

Ai 

+ A2 

+ 

Define 

ai 

= 

X3  - 

X2 

a2 

= 

Xl  - 

X3 

33 

= 

X2  - 

Xl 

bi 

= 

72  - 

73 

b2 

= 

73  - 

71 

b3 

= 

71  - 

72 

Then  the  differentiation  formulae  for  the  two-dimensional  case  are 


^ 3 3L_,  . ^ b.  „ 

q 1 d _ 3^  _i  ^ 

3x  " Zw  3x  3L  ° 2-/  2A  3L  ’ 

1=1  1 = 1 ^ 


a 3 3L,  . 3 a,. 

o _ 1 3 _i  3 

3y  2^  3y  3L.  ° 2-(  2A  3L. 

1=1  1=1  ^ 


where  Integration  over  the  area  A results  as 


p q r p!  q!  r!  2A 

L L L dA  = 

1 2 3 (p+q+r+2)! 


Note  that  the  area  of  the  triangle  simplifies  to  be 


2A  = a3b2  - a2b3  = aiba  - aabj  = a2bi  - aib2 


When  utilizing  rectangular  elements,  the  triangles  resulting  from  dividing  the 


element  have  equal  areas  of  one-half  the  rectangle's  area. 


Derivation  of  element  "XK"  and  "XC"  matrices 

The  governing  partial-differential  equation  in  the  two-dimensional  transient 
conduction  heat  problem  on  an  area  A,  with  boundary  surface  A is 


, lit  + k lit 


with  boundary  conditions  t(0)  = t^  on  the  surface  s,  a subset  of  A^;  and 


A 

s 


A1 


1 


The  initial  condition  is  t(0  = 0)  = t^  . The  derivation  of  the  two- 
dimensional  model  parallels  the  derivation  of  the  three-dimensional  model 
(see  Chapter  4),  except  that  the  z-coordlnate  term  is  omitted.  Using  the 


notation  introduced  in  Chapter  4,  let 


I = 1/2 

X 


/“x  (t) 

A 

f\  (t) 


where 


I = 1/2 

y 


V = dA 


1 = I + 1 + I 

X y pc 


For  a subarea  A of  the  original  area  A,  define 


I®  = I®  + I®  + I®  = 1/2 

X y pc 


Ait) 


+ 1/2 


r 9tj 
y p"  96  - 


I 

J 
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Then,  if  the  area  A Is  discretized  into  m elements.  Equations  11  and  12 
can  be  combined  as 


I I' 


e * 1 


m 


& e 

+ I + I 
y pc 


4M 


€.  = 1 


+ 1/2 


if  JC  ^ dA] 


a quadratic  interpolating  shape  function  is  used  to  approximate  the  field 
variable  (temperature). 

In  matrix  notation 


t = I N 1 { t^  } 


where  t is  the  temperature  (as  a function  of  position  within  an  element  "e") , 
{ t j } is  the  column  vector  of  nodal  temperatures  of  element  "e"  (in  the 
sequence  shown  in  Fig.  7),  and  [N  i is  the  row  vector  of  shape  functions 
defined  as 


Substituting  Equation  14  into  Equation  13  transforms  the  transient  heat  con- 
duction relation  within  an  element  "e"  into  a function  of  nodal  point  tempe- 
rature values.  This  process  is  repeated  for  each  element  "e"  of  the  discretized 
area  A.  The  resulting  relations  are  combined  and  minimized  with  respect  to 
each  nodal  variable,  producing  a system  of  linear  equations.  The  boundary 
conditions  (specified  temperatures)  and  initial  condition  are  inserted.  Then 
values  for  nodal  temperatures  are  computed  at  specified  time  step  intervals 
by  the  Crank-Nlcolson  time  advancement  routine.  Appendix  B contains  further 
mathematical  details. 


i 

i 
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Chapter  6 


DERIVATION  OF  PHASE  CHANGE  APPROXIMATION 


Discussion 

Phase  change  Is  the  major  heat  transport  phenomena.  A modified  form  of  the 
procedure  used  by  Bafus  and  Guymon  (1975)  Is  adopted  for  this  model. 

Consider  a nodal  point  "A"  In  a discretized  two-  or  three-dimensional 
domain.  Phase  change  Is  not  permitted  to  occur  at  node  A until  a requisite 
amount  of  latent  heat  has  been  exhausted  or  added.  The  prevailing  quantity  of 
latent  heat  for  node  A is  determined  by  considering  the  weight  of  material 
(subject  to  phase  transformation)  of  each  contributing  volume  associated  to  node 
A in  a nodal  latent  heat  accumulator  array. 

As  an  example,  consider  the  cooling  and  freezing  of  node  A.  After  each 
time  step  (or  a specified  number  of  time  steps) , the  temperature  of  node  A is 
tested  to  determine  whether  or  not  it  Is  within  a specified  tolerance  of  a 
prescribed  temperature  where  phase  transformation  Is  assumed  to  occur,  the 
freezing  point  depression,  T^  . If  the  node  has  not  been  previously  frozen 
and  the  computed  temperature  has  dropped  below  T^  , then  the  quantity  of  latent 
heat  evolved  during  the  previously  computed  time  step  Is 

G = Cu  (Tj  - COMPUTED  TEMPERATURE)  (VOLUME)  (SPECIFIC  WEIGHT) 
where  C^  « a weighted  specific  heat  (at  constant  temperature) 


A5 


Volume 


volume  of  material  assigned  to  node  A 


Specific  weight  = a weighted  average  of  specific  weights 

The  quantity  G is  subtracted  from  the  latent  heat  accumulator  for 
point  A and,  if  more  latent  heat  must  still  be  exhausted,  the  computed  temperature 
at  the  node  is  brought  back  to  prior  to  proceeding  with  the  time  advancement 
process.  This  is  analogous  to  simulating  phase  change  as  an  isothermal  process. 
The  time  rate  of  latent  heat  generation  is  governed  solely  by  the  solution  of 
the  heat  conduction  equation  and  is  generally  insensitive  to  large  time  steps, 
or  groups  of  time  steps. 

If  the  required  amount  or  more  than  the  required  amount  of  latent  heat  has 
been  extracted,  node  A and  the  volume  of  soil-water  assigned  to  it  are  considered 
to  be  frozen.  The  thermal  properties  of  the  volume  common  to  node  A are  then 
updated  based  on  the  volumetric  proportions  of  soil-water  and  soil-ice  present 
in  the  respective  elements.  If  an  excess  amount  of  latent  heat  is  removed, 
the  residual  is  used  to  calculate  how  far  below  T^  the  temperature  at  node  A 
should  be. 

The  thawing  process  is  handled  in  a similar  manner.  Provisions  are  made 
to  monitor  which  points  are  frozen  or  thawed  so  as  to  avoid  refreezlng  frozen 
areas.  Boundary  conditions  remain  unaffected  by  the  checking  routines  that 
scan  point  temeratures.  After  phase  transformation  at  a point  occurs,  the  latent 
heat  accumulator  assigned  to  that  point  is  recomputed.  Recomputation  permits 
the  simulation  of  long-term  cyclic  freezing  and  thawing. 
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Change  of  State 

Again,  consider  the  freezing  of  node  A In  the  discretized  continuum.  Once 

the  temperature  of  node  A falls  below  the  freezing  point  depression,  the  latent 

heat  accumulator  assigned  to  node  A Is  adjusted  per  the  amount  of  heat  added 

or  exhausted.  Using  the  specific  heat  parameter  for  the  unfrozen  material, 

the  amount  of  heat  evolved  In  a change  of  temperature  is 

Q = me  (t.  - t.) 

u f 1 

where  m = mass  (or  weight,  depending  on  units) 

c^  = specific  heat  (unfrozen) 


t^  = temperature  of  final  state  j 

1 

t.  = temperature  of  initial  state 
^ ( 

j 

In  the  English  system,  for  example,  the  specific  heat  of  water  (unfrozen)  j 

Is:  1 Btu/lb-*^F;  and  of  ice:  0.51  Btu/lb-°F  . The  amount  of  heat  that  must  i 

I 

be  evolved  per  unit  mass  to  freeze  (or  thaw)  water  is  approximately 

1A4  Btu/  Ibm  or  80  cal/gm  . If  more  than  the  necessary  amount  of  latent  heat  ^ 

Is  exhausted,  the  excess  heat  evolved  is  used  to  calculate  the  new  nodal  i 

temperature,  using  the  frozen-state  specific  heat  at  constant  pressure. 

The  following  example  illustrates  the  above  procedure: 

Given  that  10  lbs.  of  material  is  assigned  to  node  A.  Assume  that  this  material 
is  100%  water  and  assume  that  the  temperature  of  node  A actually  is  the  average 
temperature  of  the  material  assigned  to  it.  (See  the  next  section  for  further 
discussion  of  average  temperature).  Finally,  assume,  for  example,  that  the 
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previous  time-step  temperature  of  node  A was  , and  the  new  computed  J 

temperature  is  29°F.  | 

This  change  of  temperature  represents  a heat  loss  of  ; 

Q (Btu)  = (10  lbs.)  (1.00  Btu/lb-”F)  (29°F-34°F)  =“50  (Btu)  ‘ 

However,  node  A cannot  assume  29°F  until  the  requisite  (144  Btu/lb)  (10  lb)  = | 

1,440  Btu  of  latent  heat  is  exhausted.  Hence,  node  A would  be  set  at  32°F,  j 

and  the  latent  heat  accumulator  reduced  as  follows:  i 

LATENT  HEAT  (NODE  A)  = 1,440  + (10)  (1)  (29°F  - 32*^^  | 

=1,440-30 
= 1,410 

Consider  that  the  latent  heat  accumulator  for  node  A is  almost  zero,  perhaps 
10  Btu,  and  the  new  computed  temperature  for  node  A is,  for  example,  25°F. 

Then  we  proceed  as  follows: 

Q = (10  lbs)  (1.00  Btu/lb  - °F)  (32°F  - 25°F) 

= 70  Btu,  which  is  greater  than  the  remaining  latent  heat  for  node  A 
This  implies  the  material  assigned  to  node  A is  now  frozen,  and  yet  60  Btu  of 
heat  remains  to  be  accounted  for.  A new  temperature  has  to  be  calculated, 
where  a change  of  temperature  is  a function  of  the  specific  heat  of  ice,  at 

0.51  Btu/lb-°F,  and  the  amount  of  heat  involved.  Hence,  the  new  temperature  of  | 

node  A is  computed: 

60  Btu  = (10  lbs.)  (0.51  Btu/lb-°F)  (32°F  - t^) 

Thus,  11.76  = 32  - t^  or  t^  = (32  - 11.76)  ^F  = 20.24°F  ; 

The  latent  heat  accumulator  for  node  A is  now  zero,  and  the  time  progression  i 

■i 
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r ^ 


continues,  but  with  node  A now  set  initially  at  20.2A°F  . Note  how  the  phase 
change  dominates  the  process  of  lowering  the  temperature  of  node  A from  34°F 
to  20.2A°F. 

Average  Temperature  - Two  Dimensions 

Each  triangular  element  contains  six  nodes.  Each  node  is  a sample 
temperature  within  the  region  of  material  associated  to  it.  It  will  be 
shown  how  the  material  is  associated  to  each  node,  and  how  the  nodal 
temperature  is  used  to  determine  phase  change. 

Each  triangular  element  can  be  considered  as  the  union  of  four  triangles, 
each  containing  three  nodal  points.  The  boundaries  of  the  geometry  associated 
to  the  nodes  are  determined  by  the  perpendicular  bisectors  of  the  sides  of  the 
interior  triangles  (see  Fig.  10).  The  model  makes  the  following  assumption; 

When  a nodal  temperature  is  within  a specified  tolerance  of  the  phase  change 
temperature,  the  nodal  temperature  is  assumed  to  represent  the  average  temperature 
of  all  the  associated  material  to  the  node,  and  this  temperature  is  then  used 
for  phase  change  determination  as  discussed  in  the  previous  section.  But  is 
this  simplification  valid?  The  following  discussion  approaches  the  problem  by 
a rigorous  interpretation. 

Consider  the  freezing  of  node  A,  as  illustrated  by  the  hatched  geometry  in 
Fig.  10,  and  a computer  model  utilizing  exact  temperature.  The  computer  program 
keeps  track  of  each  element  containing  node  A.  Once  node  A is  tested  for  phase 
change,  each  element  determined  as  containing  node  A is  scanned  as  to  the 


A 
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A quadratic  element  is  consiierel 
triangles-each  with  three  nodal  p 


3ach  interior  triangle 
is  subdivided  by  the  extensions  of 
perpendicular  bisectors  of  each  si* 
of  the  triangle 


the  union  of  four  interior 


A continuum  subdivided 
the  into  the  geometries 

e associated  to  each  node 

for  phase  change  study. 


Fig,  10.  Phase  Change  Scheme 
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location  of  node  A within  the  element  ( 1.  e.  as  to  whether  node  A Is  a midside 


node  or  a corner  node).  The  area  of  each  element  associated  to  node  A by  the 
geometric  subdivision  through  perpendicular  bisection  (Fig.  10)  of  the 
triangle  sides  Is  calculated,  then  multiplied  by  the  average  temperature 
associated  to  this  sub-area,  and  then  this  product  or  weighted  contribution 
Is  added  to  the  similar  weighted  contribution  for  each  element  associated  to 
node  A. 

This  sum  of  weighted  contributions  Is  divided  by  the  total  area  associated 
to  node  A,  resulting  In  the  average  temperature  of  the  geometry  assigned  to 
node  A. 

Consider  now  the  average  temperature  of  an  entire  quadratic  triangular 
element  In  our  two-dimensional  problem.  The  average  temperature  for  the  entire 
area  Is  simply 


t 

e 


( N 1 (Ji^dA 


where  A Is  the  area  of  the  element,  [N  ] Is  the  Interpolating  shape  function 
matrix,  Is  the  vector  of  nodal  temperatures. 

From  a previous  section,  we  have 


At  = /*  (2L  - L ) <{)idA  + f (2L  - L )(|)2cii^.,+ /*  (2L^  - L ) <>3dA 

eey  11  J 2 2 J 3 3- 


+ 

A 


/4L  L 0 dA  + / 4L  L (J)  dA  + / 4L  L d)  dA 

> 2 “ J 2 3 5 J 3 16 
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A/3  (<))4  + <{)s  + ^e) 


••  t = 1/3  ($4  + 't's  + 'Pi,) 

e 

the  average  of  the  midside  nodes.  Consider  the  average  temperatures  of  each 

subtriangle  in  the  top  of  Fig.  10.  The  sum  of  the  average  temperatures, 

multiplied  by  their  respective  area,  should  ultimately  sum  up  to  be  equal  to 

4At  , 
e 

4 

i.e.,  1/4A^  A^t^  = = 1/3  (ip4  + (p5  + 4>6) 

1 = 1 

We  shall  now  derive  these  average  temperatures,  tA^  ,1  = 1,  2,  3,  4 . 

Consider  the  element  oriented  as  shown  below: 
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I 


the  local  coordinate 


= y/b 


■■■  2L  - L 

3 3 


= 2y^/b^  - y/b 


then  the  integration  for  4>3  results  In 


a/2  bx/a 


m 


N 


X + 


bS, 


2a-a) 


Ast 


J J \b^  b 

(03)  0 0 


a/2  0 


-b& 

96 


, Aa  - (1/2)  (2/2)  (b/2)  = b2/8 


^‘0’^  ^Aa  = = - IT 


By  argument  that  the  triangle  Is  arbitrarily  oriented,  the  coefficient  of  02 
must  also  be  - 1/12  . 

Consider  a new  orientation  of  the  element: 
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Th«.p  the  Integration  for  i))i  results  Ln; 

a bx/a 


'‘a.  ■ / / « < a i> 


( ^ + ^))  dydx 


„ , 2bx-ba 

2a-?. 


/ 7 


( 2b  2”^  ^ 


a 0 


= 3b?/48 


where 


A3  = (i/2)  (b/2)  = bi/A 

■•■  U.IJA.-  (^) 


So,  for  an  arbitrary  corner  subtriangle  of  a quadratic  triangular  element  we 
have  intermediate  values  for  the  average  temperature: 


t^^  - 3/12  (})i  - 1/12  4)2  - 1/12  4>3  + + 624*5  + 634*5 


A A\ 

'A  / «4. 


By  analogy  that  Aj  and  Am  are  ilso  comer  subtrlangles , 


t * -1/12  <t>i  - 1/12  4)2  + 3/12  4>3  + ai4>4  + a24>5  + a34>6 

Ai 


t.  = -1/12  4)1  + 3/12  4>2  - 1/12  4)4  + Yi4'4  + Y24)5  + Y34'6 

Am 


But  since  t = 1/3  (4>4  + 4)5  + 4)6)  we  can  deduce  that 

e 


It,  = -1/12  4)1  - 1/12  4)2  - 1/12  4)3  + ni4)4  + n24)5  + n34)6 

A2 

By  the  previous  work,  we  conclude  that  an  arbitrarily  formed  triangle  can  be 

rotated  for  midside  nodal  Integrations;  hence,  we  deduce  that  Hi  = 02  = >13  • 

‘ 6 

I Therefore,  n,  = 5/12  , i = 1,  2,  3,  (also  note  that  £ (coef f icients)4)  =1 

[ ^ 1=1  ~ 

1 for  each  subtriangle). 

Consider  the  following  orientation  of  the  quadratic  triangular  element, 
and  the  subsequent  local  coordinate  calculations: 

1 

I 


(a‘  + b")  ' (x‘  + y‘ 


1/2  (a^ 


, r 

r i 

^ bxy  “ 

-0 

fL2L3dA  = 

A3 

dA,  where  As  = -^  (by  diagram) 

all  bx/a 


/ 2bx-bi.  \ 

in  \2(.a-l)) 


f I f(- 


a/2  0 


Thus,  all  midside  nodes  opposite  the  subtrlangle  being  Investigated  have  a 
temperature  contribution  of  1/12  , and  by  previous  arguments  of  orientation. 


ai  = 62  = Y3  ; Cla  = “3  = 61  = 63  = Yi  = Y2 

Thus,  ai  = 1/12  and  02  = 5/12  . 

t - -1/12  <i>i  - 1/12  ((>2  + 3/12  (t>3  + 1/12  $4  + 5/12  4>s  + 5/12  4>(, 

Ai 

t.  = -1/12  <})i  - 1/12  4)2  - 1/12  4>3  + 5/12  4)4  + 5/12  4>5  + 5/12  4)6 

A2 

t.  = 3/12  4)1  - 1/12  4)2  - 1/12  4)3  + 5/12  4)4  + 1/12  4)5  + 5/12  4>6 

A3 

t,  = -1/12  (t)i  +3/12  4)2  - 1/12  $3  + 5/12  4)4  + 5/12  4)5  + 1/12  4>6 

A4 

4 

where  1/4  ~ ~ 

1=1  ^ 

By  geometry,  each  = A/4  , for  A = area  of  quadratic  element.  We  will 

assume  that  (l/3'^A^  Is  contributed  to  each  of  the  three  nodes  attached  to  A^  , 
thus  a contribution  of  A^/12  Is  weighted  to  each  of  the  nodes  per  each 
associated  subtriangle. 

We  thus  standardize  the  weighted  contributions  of  the  overall  quadratic 
element  per  each  of  the  six  nodes,  letting  W.C.  (1)  = weighted  contribution 

for  node  (1): 
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each  corner  node  (fii,  <})2 , <}>3 


W.C.  (1)  = A/12  t,  ,1  = 1,  2,  3 

*1 

where  is  that  subtriangle 


associated  to  node  "1" 


each  midside  node  i()4,  05.  06 


W.C.  (4)  = A/12  (t.  + t.  + t.  ) 

As  Au  A2 


W.C.  (5)  - A/12  + E^^  1 Ej_) 


W.C.  (6)  - A/12  (E.  + E,  + E.  ) 

Aj  Ao  A3 


or  equivalently, 


W.C.  (4)  - A/12  (1/12  01  + 1/12  02-  3/12  03  + 15/12  04  + 11/12  05  + 11/12  06) 


W.C.  (5)  - A/12  (-3/12  01  + 1/12  02  + 1/12  03  + 11/12  04  + 15/12  05  + 11/12  06) 


W.C.  (6)  - A/12  (1/12  01  + -3/12  02  + 1/12  03  + 11/12  04  + 11/12  05  + 15/12  06) 
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wht*  le 


6 ; 
W.C.  (1)  = A^yi2  (<!.8/12  4)4  + 48/12  4)5  + 48/  12  4)6)  ’ 

i 

1 = 1 : 

= A /'i  (4)4  + 4's  + 4>6)  = A 

e e c 

which  agrees  with  the  above  derivations  of  t^  . 

We  shall  use  the  above  derivations  to  support  the  simplification  of  letting 
the  calculated  value  of  node  "i"  be  the  average  temperature  of  the  associated 
area  to  node  "i"  . 

"Aspect  Ratio"  reasoning  dictates  that  element  diameter  ratios  should  be  as 
closed  to  1 - 1 as  possible.  Thus,  one  trianglular  element  would  be  of  the  same 
approximate  shape  and  size  as  its  neighboring  elements,  otherwise  errors  would  be 
introduced  in  the  exaggerated  element's  direction  (see  Desal  and  Abel  (1972)). 

Consider  an  arbitrary  node  "1"  and  its  associated  area  (Fig.  10).  The 
average  temperature  is  the  sum  of  all  weighted  contributions  divided  by  the 

total  contributed  area.  By  the  above,  the  temperature  at  node  "A"  approximates  ^ 

the  average  temperature  for  the  associated  area,  hence  eliminating  the  above 
calculations.  The  reasoning  follows  for  a simple  midside  nodal  calculation: 
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let  Ai  % A2  , 

(Ai  + A2) 

t,,  ,,=  A/12  (-3/12  (4.1  + 4>5  1 + 2/12  [4)7  + 4'3] 

4 * 

+ 11/12  I 42  + 4>6  + 414  + 4>8  1 + 30/12  [ 4>9  1 ) 

If  the  triangles  are  similar,  we  can  see  each  bracketed  sum  approximates  4>9 
in  the  average,  hence 


which  Implies  that  *=  <p9  . 

Thus,  the  more  exacting  computer  model  (which  employs  temperature  averaging 
of  all  contributing  elements)  can  be  simplified  by  letting  the  nodal  temperature 
represent  the  average  temperature  of  the  material  associated  to  the  node.  This 
saves  a significant  portion  of  the  computer  time  involved  in  a phase  change 
process,  and  eliminates  the  necessity  of  storing  the  nodal  locations  of  each 
associated  element. 

Computer  Simulation  of  Element  Phase  Change 

The  phase  change  process  is  modeled  ' in  the  computer  program  by  updating  the 
thermal  parameters  of  any  element  whose  nodes  have  all  experienced  the 
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transformation  of  either  freezing  or  thawing.  Noting  that  the  phase  change 
portion  of  the  program  is  restricted  to  isotropic  (thermal)  domains,  the 
thermal  conductivities  in  the  x,  y and  z directions  are  all  equal  to  the 


i 

1 

i 

I 

i 


parameter  designated  as  or  "kj"  , l.e.  the  thermal  conductivity  for 

unfrozen  and  frozen  material  respectively.  Thus,  once  a node  freezes,  each 
element  associated  with  tlie  newly  frozen  node  is  tested  to  determine  whetlier 
all  the  nodes  in  the  element  are  now  frozen;  if  so,  then  the  parameters  are 
are  updated  for  that  element,  if  not,  the  program  continues  without  an  element 
phase  change  modification. 

To  modify  the  global  matrices  for  the  changed  thermal  conductivity  parameter 
of  a newly  frozen  element,  the  difference  (k^  - k^)  is  calculated  and  then 
used  as  the  thermal  conductivity  of  the  newly  frozen  element,  and  the  element 
is  reprocessed  in  the  element  matrix  derivation  subroutine.  This  technique 
alters  the  global  matrix  system  without  revising  the  entire  global  matrix  system 
element  by  element.  For  a newly  thawed  element,  the  above  procedure  is 
paralleled,  but  with  the  element's  thermal  conductivity  parameter  adjusted  by 
the  difference  (k^  - k^)  . See  Appendix  C for  "global  matrix  phase  cliange 
updating"  derivation. 
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Chapter  7 


CRANK-N I COLSON  METHOD 


Once  the  volumetric  brick/tetrahe.dron  elements  have  been  amalgamated  into 
the  global  matrices  and  where  is  the  global  conduction  matrix 
and  is  the  global  capacitance  matrix,  the  problem  becomes  one  of  solving 
the  following  system  of  linear  first  order  differential  equations: 


dt 


GKt  + Get 


0 , 


or 


GCt  = -GKt 


where  t and  t represent  the  time  differential  of  the  field  variable  vector,  and 
the  field  variable  vector  . respect! vely. 

The  Crank-Nlcolson  method  moves  the  solution  vector  ahead  in  time  Intervals 
of  A6  by  the  relation 


,(v+l)  ^ j.(v)  M |^(v)  |.(v+l)  I ^ 


or  equivalently 


GCt 


(v+1) 


GCt^''^  +—  (cct^''^  + GCt^^'*'^^)  , 

^ 2 ^ ^ 
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Combining  terms, 


(sf  + f S!s) ■ (ffi-r  £5'^''*) 

where  for 

A0  > 0 . (le  S?  + 2=  - (l5  £?  - £5)  <» 

The  global  matrices  "gK"  and  are  modified  to  represent  Equation  1 

without  additional  storage  requirements.  By  letting  the  matrix  GC*  equal 
■fo  SS  ’ form  System  1 by  setting  ® + GC*;  GC**=  2GC*  - , 

where  in  each  new  relation,  the  new  forms  of  the  GC  and  GK  matrices  are  used 
as  they  are  produced.  Hence,  the  System  1,  modified  as  above,  reads  as 

GK*t^'^^^  = GC**t^''^  (2) 

IV/  A/ 
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It  is  this  form  of  Equation  1 that  is  directly  formed  In  the  computer. 
Before  solution  of  Equation  2,  the  boundary  conditions  (specified  temperatures) 
are  Inserted  by  the  process  described  in  Meyers  (197] );  let  t^  = h,  a sperifjt  . 
boundary  temperature,  then  for  Equation  2: 


1 

& " 

11  ^ 

1 

1 

<• 

n 

9 " o j 

1 

i 4. 

! 

! 

i 

1 ! 
- fi4 
TT* 

where  the  "$"  represent  column  terms  in  the  matrix,  and  the  "c"  represent 
column  terms  in  the  GC**  matrix. 

(1)  set  (5,5)  = (5,5)  = 1.0 

(2)  set  TT*  (I)  = n ($  - C),  I 5;  -TT*  (I)  = TT  (I) 

(3)  set  "5"  row  and  "5"  column  of  gK*  & ^**  = 0.0,  except  for 
statement  (1) . 

(v+1)  Cv) 

Now  the  system  = 59’'  JJ  ’ “here  and  TT  represent 

the  modified  matrices  of  Equation  2 (but  with  inserted  boundary  conditions)  is 
ready  for  a time-advancement  solution  scheme  with  specified  time  steps. 
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Chapter  8 


MODEL  EVALUATION 


Variables  of  Convergence 

The  computer  program  utilizes  a quadratic  shape  function  of  the  finite 
element  approximation  for  temperature  (field  variable)  in  a three-dimensional 
or  two-dimensional  continuum.  As  discussed  in  Myers  (1971),  the  solution 
oscillates  about  the  analytic  solution  before  convergence.  A balance  of  element 
size,  shape  time  step  size,  and  the  magnitude  of  the  difference  between  the 
initial  temperature  and  boundary  temperatures  is  required  for  an  optimum  solution 
to  the  problem. 

Thirty  runs  were  made  varying  individually  the  above  mentioned  parameters  on 
a control  problem.  A basic  transient  heat-conduction  problem  was  employed  and 
compared  to  the  analytic  solution  by  varying  the  (1)  size  and  shape  of  the  elements; 
(2)  the  time  step  increments  and  (3)  the  temperature  "shock"  to  the  system. 
Convergence  Analysis 

A six  foot  width  "plane  wall"  problem  was  examined  for  the  following  cases: 

(1)  Approximate  the  system  by  six  cubic  elements  of  one  foot  dimensions. 

Use  a three-dimensional  model.  Shock  the  system  by  a temperature 
difference  of  200“F.  Use  a time  step  of  0.1  hours.  The  results  were 
considered  "converged"  after  approximately  5 hours  into  the  solution. 

A two-dimensional  approximation,  with  similar  parameters,  was  con- 
ducted and  revealed  similar  answers.  Thus,  the  two-dimensional  case 
was  used  for  the  remainder  of  the  testing.  The  necessity  for  quick 
convergence  to  the  exact  solution  is  primarily  due  to  the  solution 
deviations  caused  by  phase  changes  at  Incorrect  times.  Hence,  if 
convergence  criteria  is  established  for  large  temperature  "shocks” 
to  the  system,  convergence  is  assured  for  the  anticipated  temperature 
differences  imposed  in  a more  natural  situation,  such  as  5“F,  more 


or  less.  Additionally,  sudden  "freezes,"  as  may  occur  In  chilled 
climates  with  sudden  high  velocity  winds,  will  not  destroy  the 
model.  Two  shocks  for  the  three-dimensional  system  were  used,  a 
100°F  and  200®F  temperature  change.  The  two-dimensional  model  yielded 
the  same  results  In  both  tests. 

(2)  Using  the  two-dimensional  model,  a temperature  "shock"  of  two  hundred 
degrees  Fahrenheit,  and  an  element  size  of  1 foot  dimensions  (square), 
vary  the  time  step  increments  letting  "theta"  equal  0.05,  0.10,  0.25, 

0.50,  1.00,  2.50,  5.0,  10.0  and  2A.0  hours. 

(3)  Using  the  two-dimensional  model,  a temperature  "shock"  of  two  hundred 
degrees  Fahrenheit,  and  time  steps  of  0.05,  0.1,  0.25,  0.50,  1.00, 

2.50  and  5.0  hours  vary  the  size  of  elements  (square)  at  1 foot,  1.5  feet, 
and  3.00  feet. 


Analytic  Solution 

The  problem  considered  is  actual  a one-dimensional  "plane-wall"  problem. 

The  wall  is  assumed  to  be  iron  (Isotropic  with  thermal  conductivity  of  31.4  Btu/hr 
ft-°F;  density  of  492  Ibm/ft^;  heat  capacity  of  0.137  Btu/lbm-°F. 

The  governing  differential  equation  is 


3^t 


9t 


~ 96 


where  the  boundary  conditions  are  that  t « t @ x = 0,6.  The  initial  condltlc 


is  t(x,  9“0) 


t^.  The  problem  is  normalized  using  the  variables 


3C 

6 


ft*  - 

® 36 


t-t 


i on 


a= 


pc 


Then  the  differential  equations  simplifies  to: 


; t*(O,0*)  = 0,  t*(l,9*)  = 0,  t*(x*,0)  - 1 

The  solution  is  the  power  series: 

4 1 -9tt^0*  1 -25'it^0* 

t*(X*,0*)  = — (sin  x*e  + — sin  3irx*e  + — sin  5TO*e  +. . .) 
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Solutions  to  this  equation  for  the  case  where  the  boundary  condition 
temperature  at  x = 0 is  lOO^F  and  at  x = “°  is  - 100“F  are  shown  in  Table  1. 


Test  Results 

The  model  output  was  examined  as  to  when  convergence  occurred.  The  time 
(in  model  simulation)  of  convergence  was  plotted  graphically  as  a function  of 
time-step  size  and  element  size.  For  the  Transient  Heat  Conduction  model,  a 
time  step  size  of  0.1  hours,  and  the  convenient  element  size  of  1.0  foot  was 
found  adequate  in  rate  of  convergence.  Results  are  shown  in  Figures  11,  12  and 
13. 


3-D  Test  Problem 


The  local  coordinate  relation  to  global  coordinates  is  stated  in  Desai 
and  Abel  (1972')  as 


c '\ 

i 

1 

1 

\ 

1 

X 

r- 

xC)') 

'ikCs 

xC-i) 

Y 

YCi') 

XW) 

iCt') 

L,' 

Lv 


for  the  tetrahedron 


The  corner  coordinates  are  numbered  by  choosing  any  vertex  as  (X,Y,Z), 

and  numbering  the  remaining  three  corners  in  a clock-wise  direction  as  viewed 
from  Pj^.  The  coordinates  of  the  P^  are  then  substituted  into  the  above  matrix. 

The  coordinate  system  is  established  in  the  figure.  Note  that  the  y- 
coordlnate  is  positive  into  the  paper.  Thus,  the  volume  under  consideration 
must  be  dimensional  from  the  origin  shown  above.  The  problem  tested  in  the 
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routine  is  an  insulated  iron  homogeneous  isotropic  bar.  The  bar  is  six  feet 
long,  with  a square  cross  section  of  one-ft^  area.  The  bar  Is  Initially 
assumed  to  be  at  lOO^F.  Then,  the  ends  are  both  set  at  0“F.  Due  to  symmetry 
the  problem  is  simplified  to  an  insulated  bar  of  1 ft^  cross  sectional  area, 
insulated  at  one  end  and  along  the  surface  boundary,  initially  set  at  100°F, 


with  one  end  set  at  O^F. 


(or  -\cxiPf^ 


-O'F 

Cor 


The  volume  is  discretized  into  3 cubic  brick  elements  of  equal  volume  v' 
The  elements  are  formed  below: 


ii  50  ae  38 


44-  531 


0 °46 


24  24  51 


44  49 


element  1 


element  2 


SB 

cilement  3 


Parameters  are: 

NBAND  = 26  (not  necessarily  exact,  but  must  be  > actual  bandwidth 
MODES  = 60  (number  of  nodes) 

NELE  = 3 (number  of  elements) 

NBRICK  = 3 (number  of  brick  elements) 

NTETRA  = 0 (number  of  tetrahedron  elements) 

NUMBC  = 9 (number  of  specified  temperature  nodes) 

THETA  = time  increment,  in  hours 
DAYS  = simulation  duration,  in  hours 

XKX,XKY,XKZ  = 31. A (thermal  conductivity,  in  Btu/hr.-ft.-‘’F 
(Note  that  the  brick  element  nodal  sequence  input  ia:) 

ELEMENT  1:  1,10,18,19,20,12,3,2,11,4,13,21,22,23,14,6,5,7,15,24,25,26,17, 
9,8,16. 

ELEMENT  2:  18,27,35,36,  etc. 

ELEMENT  3:  35,44,52,53,  etc. 


P 
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TABLE  1 


ANALYTICAL  SOLUTION  TO  TEST  PROBLEM 


Time 

hours 


1 

40.01 

92.29 

99.25 

2 

-1,37 

70.83 

88.74 

3 

-13.19 

50.36 

73.62 

4 

-23,60 

32.33 

52.81 

5 

-32.75 

16.47 

34,49 

6 

-40.81 

2.51 

18.37 

7 

-47.91 

-9.78 

4.18 

8 

-54.15 

-20.59 

-8.31 

9 

-59.65 

-30.11 

-19.30 

10 

-64.48 

-38.49 

-28.97 

15 

-81,24 

-67.51 

-62.49 

20 

-90.09 

-82.84 

-80.19 

24 

-94.06 

-89.70 

-88.11 

Temperatures,  °F,  at  Indicated  Intervals 


1 foot 


2 foot 


3 foot  (midpoint) 


TIME  STEP  SIZE  I 


5 10  IS  20  25  30 

CONVEROfNCE  IHOURSI 


1 Foot  Dimensions 
2-D  Problem 
200®  Shock 

FIG.  11 
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9 10  IS  30  25  30 

CONVERGENCE  IHOURSI 

1.5  Foot  Dimensions 
2-D  Problem 
200°  Shock 

FIG.  12 
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5 


10  15  20 

CONVCROCNCE  IHOUHSI 


25 


30 


3 Foot  Dimensions 
2-D  Problem 
200“  Shock 

FIG.  13 
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APPENDIX  A 


THREE  DIMENSIONAL  FINITE  ELEMENT  MODEL  DERIVATION 


GOVERNING  EQUATION 

The  governing  partial  differential  equation  for  the  transient  heat  conduction 
model  In  three  dimensions  is 


k 


a^t 

X Bx^ 


(1) 


where  (a)  kx,  ky,  kz  are  constant  dlrectinal  thermal  conductivities. 

(b)  0 Is  the  time  variable 

(c)  p is  constant  density 

(d)  c Is  the  constant  pressure  specific  heat 

(e)  an  initial  condition  is  given 

(f)  Boundary  values  are  specified 

(g)  at/an  = 0 along  the  surface  of  the  volume  being  studied 
(i.e.,  the  volume  is  thermally  insulated). 

Although  Equation  1 assumes  constant  physical  and  thermal  properties  through- 
out the  control  volume,  variable  properties  can  be  handled  by  the  numerical 
methods  proposed  (see  section  entitled  "Derivation  of  System  Matrices,") 
contained  in  this  appendix) . 

DISCRETIZATION 

Given  a volume  V,  discretize  V into  a finite  union  of  "m"  tetrahedra  shaped 

(e) 

elements  V where 


V «=  U,  V 
e=l 


(e) 


and  the  intersection  of  any  two  V is  an  entire  common  face,  edge  or  vertex. 
For  each  tetrahedra,  specific  nodal  points  at  each  vertex  and  at  the  midpoint 
of  each  edge.  This  will  result  in  10  nodes  for  each  element  V . Due  to  the 
discretization  of  a continuous  volume,  nearly  all  of  the  nodes  will  be  shared 
by  other  elements.  A special  requirement  of  the  discretization  process  is  that 
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two  contintiguous  elements  must  share  an  entire  face  or  an  entire  edge  rather 
than  portions  or  segments.  This  must  result  in  6 common  nodes  (face)  or  3 
common  nodes  (edge)  for  contiguous  elements.  Only  when  elements  touch  at  a 
vertex  will  there  be  just  one  common  node. 

ELEMENT  AND  NODAL  NUMBER 

Number  the  elements  and  nodes  of  the  descretlzed  volume.  Do  not  renumber 
shared  nodes  of  contiguous  elements.  That  is,  if  the  system  is  discretized  into 
two  tetrahedron  contiguous  along  a face,  there  will  be  2 elements  and  14  nodes 
in  the  model,  not  2 elements  and  20  nodes. 


DERIVATION  OF  SYSTEM  MATRICES 

It  was  shown  by  Desai  and  Abel  (1972)  that  solving  Equation  1 along  with 
its  specifications  is  equivalent  to  minimizing  with  respect  to  the  variable 
"t"  (temperature)  the  variational  statement 


where  we  assume  9t/9n  = 0 along  the  volume's  surface. 

If  we  could  minimize  Equation  2 with  respect  to  "t",  we  would  arrive  at 
an  expression  for  temperature  as  a function  of  position  and  time.  Instead  of 
attempting  to  solve  the  above  analytically,  approximations  will  be  used  which 
when  substituted  into  Equation  2 result  in  a system  of  linear  equations  with 
values  of  temperatures  as  the  unknowns. 

The  first  step  is  to  formulate  a function  for  temperature  within  the 
volume  V.  This  will  be  done  by  specifying  functions  for  temperature  within 

,(e) 


each  element  V 


These  temperature  functions  will  be  continuous  along  the 


element  Interfaces  such  that  Equation  2 will  be  defined. 

C 6 ) 

Consider  an  element  V along  with  its  ten  specified  nodal  points.  Let 


(2) 


f(ti,  t^.-.-.t^Q,  x,y,z) 
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where  "t"  is  the  temperature  at  point  "P"  having  coordinates  (x,y,z)  of  the 

volume  V,  and  "t"  is  a function  of  the  temperature  of  the  element's  nodal  points 

where  P is  contained  in  V , ( (t  , t2 , . . . , t^^p)  ) are  the  element's  nodal  temperatures 

in  the  order  shown  in  Figure  5.  This  temperature  function  will  equal  the  nodal 

temperature  at  each  nodal  point,  and  will  interpolate  between  these  nodal  points 

Cs ) 

to  evaluate  temperatures  elsewhere  within  V . 

Using  the  "local  coordinate"  system  defined  in  Chapter  4,  we  can  write  such 
an  interpolating  temperature  function  as 


t = (2lJ  - L^)t^4  (2L^  - L2)t2  + (2L^  - L3)t2 
+ (2L^  - L^)t^  + 4L  L2t3  + Au^L^t^  + 


+ 4L2L3tg  + ALgL^tg  + 4L2L^tj^Q 


or  in  matrix  notation 

t = [(2L^-L3),  (2L2-L2),  2L^-L3),  (2L^-L^),  4L3L2,4L3L3,  ^L2L3,  4L3L^, 


where  t 


(e) 


4L  L ] {t^®h 
2 ^ 


,(e) 


Equation  3 is  valid  for  each  element  V within  the  volume  V.  Also,  this 
equation  represents  a quadratic  variation  of  "t"  (temperature)  throughout  each 


-(e) 


(3) 


, providing  increased  accuracy  over  a linear  interpolating  temperature  function. 


The  second  step  is  to  evaluate  the  integral  expression  shown  in  Equation  2. 
Assuming  that 


V * U,  V 

e =1 


(e) 


and  assuming  Equation  2 can  be  rewritten  as 


,9t. 


/^t^ 


9t' 


I " ^ \ (ft) 

U\/^ 


(4) 


e.*t 
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then  Equations  2 and  4 can  be  combined  and  expanded  Into 


I ^ 2,  e)  lU , (k  (|^)  + k (~)  + k (-|i)  + pc  ) dV 

e=l  2 „Ce)  X 9y  y oy  z dz'  36 


(e  ) 

Thus  we  can  consider  each  element  V individually  in  evaluating  5,  such  that 

m , , 

I = 2,  1^°^  ( 

e=l 


where 


\ fff  (k  (|^)  + k (— ) + k (|^)  + pc  )dV 

2 . V X dx  y oy  z az  30 


The  third  step  is  to  note  that  since  we  have  an  interpolating  temperature 

(e  ) 

function  Chat  is  valid  for  each  element  V , we  can  approximate  the  partial 
derivatives  in  Equation  6 by  the  differential  formulae  established  in  Chapter  4. 
Furthermore,  these  expressions  will  be  valid  for  each  element. 

(g  ) 

To  simplify  calculations,  rewrite  the  expression  for  I in  Equation  6 as 


(e)  1 


^ k (^)  dv  + ^ ;/;  k (^)  dv  + ^ ///  k (-^)  dv  + ^ ///  pc^  dv  (?) 

\(.)  “ >■  \(e)  ' 


or  simply 


I (e)  ^ I (e)  j (e)_|_  ^ (^e) 
X y z pc 


where  the  subscripts  refer  to  the  identifying  governing  paranieter  for  each 
expression  in  Equation  7.  The  reason  for  this  is 

similar,  except  for  some  multiplying  constants  of  position  and  thermal  parameters. 

( G ) 

Hence,  we  only  need  to  find  a useable  expression  for  I , and  then  expand  the 

(g)  (g)  (g^ 

results  to  include  I and  I .1  must  be  solved  separately  due  to 

y z pc 

the  time  derivative. 
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Our  problem  is  to  minimize  Equation2  with  respect  to  "t"  (temperature). 

This  is  equivalent  to  minimizing  Equation  4 with  respect  to  temperature.  But 

our  interpolating  temperature  function  is  a function  of  nodal  temperatures. 

Thus  we  must  minimize  Equation  5 with  respect  to  each  nodal  temperature  of 

( 0 \ 

the  entire  volume  V.  By  Equation  6,  we  see  that  we  must  minimize  each  I 

with  respect  to  each  nodal  temperature  of  V.  Since  only  ten  nodes  are 
(e) 

contained  in  each  V , most  nodal  temperature  minimization  efforts  (for  a 
several-element  nodel)  equate  to  zero.  That  Is,  minimizing  Equation  4 with 


respect  to  each  of  the  volume's  V nodal  temperatures  is  equivalent  to  the 
of  t 

,(e) 


sum  of  the  minimizations  of  I with  respect  to  each  nodal  temperature 


of  V' 


(0  ) 

Hence,  for  each  V we  want  to  evaluate 


9t 


(e) 

- = 0, 


(e) 


3t, 


= 0,..., 


31 


(e) 


3t 


= 0 


10 


where  t^  is  the  nodal  temperature  as  expressed  in  Equation  3. 

Since 

l(e)  ^ j.(e)  ^ j.(e)  ^ j.(e)  ^ ^(e) 

X y z pc 

(g)  (e) 

we  can  minimize  I with  respect  to  (t  , t then  extend 

X X XU 

these  results  to  the  minimization  of  I and  I . However,  I must  be 

y z pc 

minimized  with  respect  to  (t^, t2, . . . , t^^)  of  V separately. 

Step  four  is  performing  the  above  described  operations.  We  will  use  the 
notations  a^,b^  and  c^  for  the  cofactors  of  the  element's  coordinate  matrices 
(see  Chapter  4)  and  Equation  3 for  the  interpolating  temperature  function. 


Consider,  for  any  V 


(e) 
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= 


■Z  Wx  III  tiV 

“ M [c>,LUfcL^  -8L,-^0V4(A\.,U-l>s 
+ 4(A-L,Li“L^'t4  4(4.Li\-4-L4.')t,7l 

+■  a^.  LCiGL.Lz,-  4L -4Lj.-^\')t2.+  4(4 

+ 4 C4L,L3-l.3)'t8 -V"  4C4L4\-4,-L/^'t,)ol 
u.3[(\Gl..U-4LrA-L3-"'')^3+4C4Lt-0u 
^4(A.L,L^-0■t8■^  4-C4LjLa.- 14^19! 

Ol4^\GL|\_4.- 4L^“4L4+\'H4.+  4(4^“\-i')'t7 
■v-4C4LvLi'L4)‘t9  4C4L^L^-L^)■b^ol]ci^/ 


where  the  Integration  is  solved  as  follows: 


m Lt  dN  - 
it')  lJ|^L;dV= 


ui)  Snulj<i\/=(<i>^'Xito') 

iOBUv  = 


where, 


+ a^^'5b^^+  -fots  'iots  + 
+ a4(j^-i4-^'i'^'7‘^  '30^'^ 
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The  other  partial  jerivatives  follow.- 


1 

! 

5tj.  ' (4l^-l')dV 

f V‘  X \ "*  ^ “*  1 

1 c» ,L io"^'  ■*“  io^fe ■'■ 

50^8+ 

30^ 30^'^! 


•^3  ^ce) 

^ {.^'[.36'^'  ^ ^ -^io^  30"^ ■’1 

30'^®'^ 

■'■  30^"^  "io"^^  36**^‘^li 
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The  preceedlng  processes  of  partial  differentiation  must  be  repeated  for 


I 


(e) 


and  I . However,  parallel  development  will  result  in  the  same 
y z 

(e)  (e) 

expressions  as  derived  for  I , except  that  for  I we  must  substitute 


"b."  for  "a,"  , and  "k  " for  k 5 and  for  I 
i 1 y X z 

for  "a,",  and  "k  " for  "k 
i z X 


(e) 


we  must  substitute 


(e)  T 

Let  - the  column  vector  [ t t^,  . . . , tj^^]  . Then  minimizing 

I + I + I with  respect  to  (t, , t-,...,t  ) of  is  equivalent 

X V z 1 z iu 


to 


, . I + I + I 

9.C-J5 Y z 


3t 


(e) 


= 0 


(9) 


By  the  above  derivations,  we  se<_  that  Equation  9 establishes  the  matrix 
expression 

[X  K]  ^ 0 

( 6 ) (s ) 

where  [X  K]  is  the  "element  conduction  matrix"  for  V where  t'  ' is  the 

C G ) 

column  vector  of  nodal  temperatures  of  V . 

Note  that  we  can  expand  Equation  10  as 

(X  K]  = {[  = 0 


(10) 


where  each  matrix  in  the  above  sum  corresponds  to  the  minimization  of  its 

respective  entry  in  Equation  9.  The  matrix  [ is  written  in  modified 

^0  ^ C 0 ^ 

form  in  Figure  A-1  of  the  appendix.  The  matrices  [ K ] and  [ K ] can 

y z 

be  found  by  substituting  the  constants  b.  for  the  a.  and  k for  k (for  matrix 

1 i y X 

( K and  by  substituting  c.  for  the  a.  and  k for  k (for  matrix  [ K 1 . 

y 1 1 z X z 

Equation  10  only  represents  part  of  the  expression  for  minimizing 
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■•Equation  7.  V.'e  inu'^t  still  min  ir.i  z^-  with  respect  to  C'tijti,  — ,t/o'^ 

of  . that  is  we  must  solve 
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Cr:  (where  by  symmetry  f^\-  ‘=^4') 


5 


o(, 

o<3 

•^3 

•^3 

ft'A 

*^3 

0(4. 

a<4 

0^3 

oi3 

o<2. 

o(l 

■A- 

•^3 

*^4 

•^3 

•<3 

•<4 

•<4 

=^3 

oi^ 

•^3 

<^3 

A. 

ft 

^4 

1^3 

^4 

ft 

(V 

As 

ft 

(^4 

ft 

1^3 

ft 

|8, 

^3. 

(^4 

1^6 

(^3 

A4 

A 

ft 

ft 

ft 

A3 

P>4 

)3. 

^4. 

ft 

^4 

I^A- 

ft 

As 

respect 

: to 

(-b 

, -tie 

of  N 

(“  : 

• (£) 
t 


(e^ 


in  the  "element  capacitance  matrix"  expression 


(11) 


What  we  wanted  to  do  was  minimize  aquation  ? with  respect  to 
(■tvjti,  ,-tio^  Of  V Hence,  combining  aquations  ?,  10,  and  11 

we  get 

[x  ^ (12) 

By  Equations  5 and  6,  we  must  perform  the  above  operations  for 
each  of  the  "m"  elements,  and  sum  the  results  in  matrix  form.  This 
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combination  of  element  matrices  will  form  the  "global  matrix"  system 


[G  K]t  + [G  C]t  = 0 (13) 

where  [G  K]  = the  global  conduction  matrix 

[GC]  = the  global  capacitance  matrix 

t = the  column  vector  of  all  nodal  temperatures 

t = the  column  vector  of  time  derivatives  of  all  nodal 

temperatures 

Using  Equation  13,  we  can  apply  a time  advancement  modification  such  that 

the  "Crank-Nicolson"  method,  and  formulate  a matrix  system  that  can  be 

programmed  into  the  computer.  From  the  above,  it  can  be  seen  that  variable 

physical  and  thermal  parameters  can  be  handled  by  specifying  constant 

(e ) 

parameters  per  each  element  V . 

I 

1 
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r 


\8D\/“^ 


3af 

-aAt 

-a.aa 

-a^A 

(3a, 4z.-  a.^,) 

-atO., 

3 Oil 

-aj.aj 

-a.^.a4 

C3a^a,  - aO 

“O-iQ-z 

3.a| 

— O-jCL^ 

(—a^Qi,—  oljClI') 

-Q4Q.1 

-0-4Q-3 

3a  1 

(-a^a,-a4.0 

(3ai.a,'aO  C3.ajai-a.O  (-aA^-a^ci^  (-a,a^-ai,a4^(8.Q^-*-8a.at  + SQx) 


(SajCLi-OL^^  I 2)Q.Q.j- O.j')  (-Q,j(Xv“Q.^O.^XA-(i^  + 4-<X.,Cki 

44<X4a.  + 8aA.a2^ 

(■5cL4.ar  at'^  C-a.a^-a^.ai')  (-a,Ql3-aiva5^  (sa.o^-  (4^^  + A-aAi. 

4 4a4a,  vB(X4.ai) 

*^3^M3>(X3Q-i- Ck-z.')  (^•XzCk.a-'  O-j')  C~  <Xi.O4.-0ij0l^  (^djAi  4-  4-CLi 

+8aja,  + 4aj,at^ 

CL^')  (soijQ.t- Qi^')  C4a.^a.i  4 4(Z,jai 

■^Aa^CiA.-^A-a^at.) 
C4-aia, -v 4ai; 

= 4e,Q^a,4-  A-a4ai') 


.X, 


FIGURE  A-1  (appendix  A) 


92 


COLUMU'b 
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appendix  b 


TWO  DIMENSIONAL  FINITE  ELEMENT  MODEL  DERIVATION 

GOVERNING  EQUATION 

The  governing  partial  differential  equation  for  the  transient  heat 
conduction  model  in  two  dimensions  is 

3^t 


Id. 


+ k 


9t 


y 9y 


(1) 


where  (a) 

k are  constant  directional  thermal  conductivities 

y 

(b) 

0 

is  the  time  variable 

(c) 

P 

is  constant  density 

(d) 

c 

is  the  constant  pressure  specific  heat 

(e) 

an  initial  condition  is  given 

(f) 

boundary  values  are  specified 

(g) 

3t/3n  * 0 along  the  boundary  of  the  area  being 

(l.e.,  the  subject  area  is  thermally  Insulated) 

Although  Equation  1 assumes  constant  physical  and  thermal  properties  through- 
out the  control  area,  variable  properties  can  be  handled  by  numerical  methods 
proposed  (see  section  entitled  "Derivation  of  System  Matrices",  conta^’'“^d  in 
this  appendix) . 

DISCRETIZATION 


elements  A where 


A * U A' 
e*l 


Given  an  area  A discretize  A into  a finite  union  of  "m"  triangular  shaped 

(e) 

and  the  intersection  of  any  two  A'  ' is  an  entire  common  edge  of  the  triangles 
or  a common  vertex.  For  each  triangle,  specify  nodal  points  at  each  vertex 
and  at  the  midpoint  of  each  edge.  This  will  result 
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In  6 nodes  for  each  element  A . Due  to  the  discretization  of  a continuous 
area,  nearly  all  of  the  nodes  will  be  shared  by  other  elements.  A special 
requirement  of  the  discretization  process  Is  that  two  contiguous  elements 
must  share  an  entire  common  edge  rather  than  a portion  of  the  line  segment. 

This  must  result  In  3 common  nodes  for  contiguous  elements.  Only  when  elements 
touch  at  a vertex  will  there  be  Just  one  common  node. 

ELEMENT  AND  NODAL  NUMBERING 

Number  Che  elements  and  nodes  of  the  discretized  area.  Do  not  renumber 
shared  nodes  of  contiguous  elements.  That  Is,  If  the  system  Is  discretized 
Into  two  triangles  contiguous  along  an  edge,  thre  will  be  2 elements  and  9 
nodes  In  the  model,  not  2 elements  and  12  nodes. 

DERIVATION  OF  SYSTEM  MATRICES 

It  was  shown  by  Desal  and  Abel  (1972)  that  solving  Equation  1 along  with  Its 
specifications  Is  equivalent  to  minimizing  with  respect  to  the  variable  "t" 
(temperature)  the  variational  statement 

'‘y  ||’>  ■**  <2) 

where  we  assume  dt/dn  = 0 along  the  area's  boundary. 

If  we  could  minimize  Equation  2 with  respect  to  "t,"  we  would  arrive  at 
an  expression  for  temperature  as  a function  of  position  and  time.  Instead  of 
attempting  to  solve  the  above  equation  analytically,  approximations  will  be  used 
which  when  substituted  Into  Equation  2 result  In  a system  of  linear  equations  with 
values  of  temperature  as  the  unknowns. 

The  first  step  Is  to  formulate  a function  for  temperature  within  the  area  A. 
This  will  be  done  by  specifying  functions  for  temperature  within  each  element 
A'  . These  temperature  functions  will  be  continuous  along  the  element  edges 
such  that  Equation  2 will  be  defined. 


J 


1 
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^6  ) 

Consider  an  element  A along  with  its  six  specified  nodal  points.  Let 
t = f(tj,t2,...,t^,x,y) 

where  "t"  is  the  temperature  at  point  "P,"  having  coordinates  (x,y),  of  the  area 

A,  and  "t"  is  a function  of  the  temperature  of  the  element's  nodal  points,  where 

"P"  is  contained  in  A , ((t,,  t„,...,t,)  are  the  element's  nodal  temperatures 

1 Z 0 

in  the  order  shown  in  Figure  7) . 

Using  the  "local  cordinate"  system  defined  in  Chapter  Five,  we  can  write  such 
an  interpolating  temperature  function  as 

t = [(2L2-L^)t^  + (2L|-L2)t2  + (2L2-L2)t2  + ^ ^^34*^6^ 

or  in  matrix  notation  t = [ N ]t  , or 

t = [ (2L^^-L^),  (21^2-4^.  (2L3-V  >^^2’  ^4^1^ 

where  is  the  column  vector  of  the  nodal  temperatures  of  A^®\  This 

temperature  function  will  equal  the  nodal  temperature  at  each  nodal  point,  and  will 
interpolate  between  these  nodal  points  to  evaluate  temperatures  elsewhere  within 
A<^>. 

(g) 

Equation  3 is  valid  for  each  element  A within  the  area  A.  Also  this 

equation  represents  a quadratic  variation  of  "t"  (temperature)  throughout  each 
(e) 

A , providing  Increased  accuracy  over  a linear  interpolating  temperature 
function. 

The  second  step  is  to  evaluate  the  integral  expression  shown  in  Equation  2. 
Assuming  that 

A = U A^®) 
e=l 

and  assuming  Equation  2 can  be  rewritten  as 


dc 


(ft) 

;a<«> 

e“l 


(4) 


then  Equations  2 and  4 can  be  combined  and  expanded  Into 


1 " »x€>  *‘='ir>  “* 

(s  ) 

Thus  we  can  consider  each  element  A Individually  In  evaluating  5,  such  that 


(5) 


I ^ S I 

e=l 


(e) 


(6) 


where 


The  third  step  Is  to  note  that  since  we  have  an  Interpolating  temperature 

^ G ^ 

function  that  Is  valid  for  each  element  A , we  can  approximate  the  partial 

derivatives  In  Equation  6 by  the  differential  formulae  established  In  Chapter  5. 

Furthermore,  these  expressions  will  be  valid  for  each  element. 

(c) 

To  simplify  calculations,  rewrite  the  expression  for  I In  Equation  6 as 


(e)  . 1 


,3t 


r9t 


9t" 


2 ^ I V..  vw  ^ 1 fr 


,(e) 


,(e) 


(e) 


or  simply 


l(e)  _ i(e)  I ^ j 

X y pc 


(8) 


where  the  subscripts  refer  to  the  Identifying  governing  parameter  for  each 

Ce)  fe) 

expression  In  Equation  7.  The  reason  for  this  Is  I and  I ' are  similar, 

X y 

except  for  some  multiplying  constants  of  position  and  thermal  parameters.  Hence, 

(c) 

we  only  need  to  find  an  useable  expression  for  1^  , and  then  expand  the  results 

Cg)  Cg) 

to  Include  I . I must  be  solved  separately  due  to  the  time  derivative, 
y pc 
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Our  problem  la  to  minimize  Equation  2 with  respect  to  "t"  (temperature). 
This  is  equivalent  to  minimizing  Equation  4 with  respect  to  temperature.  But 


our  interpolating  temperature  function  is  a function  of  nodal  temperatures. 

Thus  we  must  minimize  equation  5 with  respect  to  eacn  nodal  tempoerature  of  the 

C 0 ) 

entire  area  A.  By  Equation  6,  we  see  that  we  must  minimize  each  I with 
respect  to  each  nodal  temperature  of  A.  Since  only  six  nodes  are  contained  in 
each  A , most  nodal  temperature  minimization  efforts  (for  a several- element 
model)  equate  to  zero.  That  is,  minimizing  Equation  4 with  respect  to  each  of 

the  area's  A nodal  temperatures  is  equivalent  to  the  sum  of  the  minimizations 

(e)  (e) 

of  I with  respect  to  each  nodal  temperature  of  A'  . 


(e) 

Hence,  for  each  A we  want  to  evaluate 
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(e) 


8t, 


= 0, 
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(e) 


91 


(e) 


9t. 


0 9t,  = 0 

O 


where  t^  is  the  nodal  temperature  as  expressed  in  Equation  3. 


Since 


l(e)  ^ j(e)  ^ j(e)  ^ ^(e) 
X y pc 


(9) 


Cs ) C0 ) 

we  can  minimize  with  respect  to  (tj^,t2, . . .,t^)  of  A , and  then  extend 

( S ) (© ) 

these  results  to  the  minimization  of  I . However,  1 must  be  minimized 

y pc 

) 

with  respect  to  (tj^,  t^)  of  A separately. 

Step  four  is  performing  the  above  described  operations.  We  will  use  the 
notation  a^  and  for  the  cofactors  of  the  element's  coordinate  matrices  (see 
below)  and  Equation  3 for  the  interpolating  temperaturefunctlon. 

Definition: 

Let 


®1  “ *3"^2 


^2-^3 


X1-X3 


b2  “ y3-yi 


83  = X2  - x^ 


yi-y2 
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where  (x^,  y^)  are  the  coordinates  of  nodal  point  location  "i”  of  the  two 


dimensional  triangular  element. 


Consider  for  any  ^ 


3t,^  dA 


then 


= C)lv^OvUv.Xr 

+ '^lOfeULr^'Lc  4.L,+iWM4Lf-LX+A-U\.,tj-\.^t^' 
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where  integration  i : nolvei  an  f'ollov/3: 


Hence , 


The  other  oartial  derivatives  follow: 


^ w r 

^ WV t5+|ttS\ 


3t-tA. 


V.K 


{ 3 ■^iC'tolvV,,V,j^v\5^")-t^ 

k C'o^^  'o,  2.v>,^  ^ 'oiW  ( v>j 

3 C\o5^iV2.\3|-'-  2.t.\  ■'’\3jSD3^tii  •'’  ^C'b3\ii+  ^l-^2.V>i'V>i+V3xV>3)ife^ 
rbati 


■^■k  C2^o.\o^+\^>3+V>sV>^  + 
3C2\oi-v\d,\>3-^ 
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r 


The  preceedlng  processes  of  partial  differentiation  must  be  repeated  for 

I . However,  parallel  development  will  result  in  the  same  expressions  as 
(6)  (c) 

derived  for  I , except  that  for  I we  must  substitute  a.  for  b, , and  k 
X y 1 1 y 


for  k . 

X 


Let  the  column  vector  [t^^,  t^,  • • • Then  minimising  + ^y^^^ 

(e) 

with  respect  to  (t, , t„ , . . . , t,)  of  A Is  equivalent  to 
Is  D 

3(1  + I 

X y r. 


By  the  above  derivations,  we  see  that  Equation  9 establishes  the  matrix 


expression 


[X  t^®^ 


where  [X  ] is  the  "element  conduction  matrix"  for  where  is 

(e) 

the  column  vector  of  nodal  temperatures  of  A 
Note  that  we  can  expand  Equation  10  as 

[X  = {[  K + [ K = 0 

X y 

where  each  matrix  in  the  above  sum  corresponds  to  the  minimization  of  its  respective 

(e' 

entry  in  Equation  9.  The  matrix  [ ' ] is  written  in  matrix  form  in  the  above 

(c) 

derivation.  The  matrix  [ ] can  be  found  by  substituting  the  constants 

a.  for  the  b, , and  k for  k . 
i 1 y X 

Equation  11  only  represents  part  of  the  expression  for  minimizing  Equation  7. 

We  must  still  minimize  l^^'^  with  respect  to  (tj^,  t2, . • • . t^)  of  A , that  is  we 
must  solve 
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Ccn^:i1er  for  any  = T ^ j®=-t\(r.e,i 

"S  >f'^elLt^f\^L^L:)c\^ 

The  other  ^partial  derivatives  J'ollow: 

^ ‘ 2(pc^'f''s  * 'Sb'ts* 


3lS 


¥ 


Ct> 


c^ 


9^^ 


“ zCpc  ^ hS>  3 


tt^ 


or  in  matrix  notation 


6.0 

-J_ 

560 

360 
O 
-i_ 

^o 

o ^ 


360 

X 

60 

-\ 

360 


:i 

360 

. 

360 


90 


O 

90 


60 

td. 

90 

O 


o 

rL 

90 

jfe 

h 

h 


o 

o 


h 

A- 

4^ 

J:. 

45 


t 

4? 

4 

2C5 


-t^ 


4i.'i 
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b(V}i+\ovV>x'''^') 

4-(t^  +bib3 

A^b.biV  ib^bj 

■^b,bi.+'2b,b^ 

+bt^'b,b5^ 

2.\5,b3'^  bib^ 

b(b|+  bib3 

^-('o.bi'Vbi 

+b.W^bl'> 

^bl^ 

■*■  Z-bbi+bib:^ 

Ubf-v  b,bt. 

4-(  b|  ♦ 2b,bi 

8(b’;+b,b%+bf) 

^ b,b3+  b^b^ 

■'•b.bs'^bit^ 

Minimizing  I with  respect  to  (t, , t~,...,t,)  of  results  in 

DC  1 ^ O 


the  "element  capacitance  matrix"  expression 


[XC^®^]  = 0 


Wliat  we  wanted  to  do  was  minimize  Equation  7 with  respect  to  (t^^, 


of  a''  . Hence,  combining  Equations  7,  10  and  11  we  get 


[X  + [XC^®^]t^®^  = 0 


By  equations  5 and  6,  we  must  perform  the  above  operations  for  each  of  the  *m' 


elements,  and  sum  the  results  in  matrix  form.  This  combination  of  element 


matrices  will  form  the  "global  matrix"  system 


[G  K]£  + [GC]^  = 0 


where 


[G  K] 


the  global  conduction  matrix 


= the  global  capacitance  matrix 


= the  column  vector  of  all  nodal  temperatures 


the  column  vector  of  time  derivatives  of  all 
nodal  temperatures 


Using  Equation  13,  we  can  apply  a time  advancement  modification  such  as  the 


"Crank-Nicolson"  method,  and  formulate  a matrix  system  that  can  be  programmed 


into  the  computer.  From  the  above  it  can  be  seen  that  variable  physical  and 


thermal  parameters  can  be  handled  by  specifying  constant  parameters  per  each 
element  A 


APPENDIX  C 


ELEMENT  PHASE  CHANGE,  GLOBAL  MATRIX  UPDATING 

When  is  the  element  considered  changed  of  phase?  In  the  proposed  model, 
nodes  are  allowed  to  change  chase  Independent  of  the  associated  element's 
phases.  Additionally,  we  can  have  frozen  nodes  in  a thawed  element,  or  thawed 
nodes  in  a frozen  element. 

The  node  changes  phase  when  the  requisite  amount  of  latent  heat  is 
appropriately  supplied  or  evolved.  However,  such  a scheme  for  an  element 
phase  change  process  is  Inadequate  in  the  critical  region  of  the  freezing 
front  where  all  nodes  could  be  considered  "frozen,"  and  yet  actually  be  in  the 
transitory  stage  of  phase  change  where  less  than  the  total  amount  of  latent  heat 
necessary  to  be  completely  frozen  has  been  evolved,  which  could  keep  the  element 
thawed.  The  model  uses  the  simple  scheme  of  keeping  track  of  the  number  of 
nodes  thawed  for  each  element. 

Each  element  is  flagged  to  indicate  whether  the  element  is  thawed  or 
frozen.  Additionally,  a vector  NTHAWD  is  used  to  store  the  number  of  nodes 
in  the  thawed  state  for  each  element. 

On  input,  the  phase  of  each  element  is  flagged  on  a stored  thermal 
parameter.  For  tetrahedron  or  triangular  elements  the  latent  heat  per  unit 
mass  is  flagged;  for  Brick  elements,  the  NTHAWD(i)  vector  is  flagged 

(a)  a positive  parameter  Implies  the  element  is  thawed  (contains 
a positive  amount  of  latent  heat) 

(b)  a negative  parameter  Implies  the  element  is  frozen. 

The  Crank-Nicolson  method  modified  the  global  matrix  system.  To  update 
this  modified  system,  the  phase-change  element  matrices  must  be  likewise 
modified.  Using  previous  notation,  the  global  conduction  matrix  G^K  and  the 
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global  capacitance  matrix  are  modified  as 


GS*  - GK  + I5  QC  I 

GC*  - I5  « - « I 

I 

t 

where  after  computation,  GK*  and  GC*  replaces  the  matrices  GK  and  QC  in 
computer  memory. 

By  definition,  the  global  matrices  C^C  and  GC  are  sums  of  the  appropriate 
element  matrices  and  The  following  example  demonstrates  how 

to  update  the  global  matrix  system  due  to  a phase  change  of  an  element,  without 
rederlving  the  entire  system  matrices.  Let 

QK  = A + 6 + £ 

OP  =£  + £ + £ 

Suppose  that  during  the  simulation,  an  element  changes  phase,  transforming 
element  matrices  £ and  £ into  matrices  £ and  Y respectively.  The  global  matrix 
modifications  from  implementation  of  the  Crank-Nlcolson  method  would  have 
previously  changed  GR  and  GC  such  that 

QIC*  = (A  +B  +C)  + Iq  ( E + £ + F) 

GC*  - (-A  -£  -C)  + Iq  (e  + £ + £) 

where  G£*  and  Q£* replaced  G{C  and  in  computer  memory.  Duplicate  the  above 
transformations  on  the  element  matrix  changes 

(X  - £)*  - (X  - £)  + le  (X  - V 
a - £)*  = (-X  + £)  + le  (X  - X) 

Now  add  (X  - C)*  and  (£  - £)*  to  GK*  and  QC*respectively  to  the  above 
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GK*  + (2  - C)*  - Oft  + £ + X)  + ^ (B  + £ + X) 

Qj:*  + (X  - £)*  = (-X  - - X)  + le  (C  + £ + X) 

which  shows  that  the  global  matrix  system  has  been  properly  modified  to  accommodate 
the  updating  of  an  element  matrix  due  to  phase  change. 

The  following  example  of  a two  element  system  undergoing  "phase  change"  will 
demonstrate  the  above. 

Let  the  modified  element  matrices  for  the  thawed  phase  be 


'l  2' 

r-A  3] 

XK*  = 

xc*  = 

.2  1. 

[3  -aJ 

and  for  the  frozen 

phase 

[3  Al 

00 

CM 

I 

XK*  = 

xc*  = 

rr 

u> 

1 

.8  2] 

Assume  both  elements  are  initially  thawed,  then  the  system  Is 


i+1 


12  0 
2 2 2 
0 2 1 


-A  3 0 
3-8  3 

0 3 -A 


where  "1"  refers  to  a time  step  Interval,  and  t^,  t^  are  nodal  temperatures. 
Assume  that  a boundary  condition  is  t^  = n-  Then 


1 

0 

0 

t 

1 

1+1 

’1 

0 

o' 

» • 

n 

0 

0 

2 

2 

^2 

= 

0 

-8 

3 

^2 

+ 

-n(2-3) 

0 

2 

1 

0 

3 

-A 

a 

^2. 

-n(o-o) 

(1) 


GK* 


GC* 


TT* 
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Now,  consider  what  the  system  would  be  had  element  number  one  been  Initially 
frozen.  The  system  would  define  the  following  relation: 


i+1  1 


— • 

o 

1 

'2  8 o' 

‘4' 

4 4 2 

*^2 

m 

8-2  3 

^2 

CM 

O 

■*^3' 

_0  3 -4. 

-h. 

where  “ h 


o 

o 

in. 

■^l‘ 

i+1 

'l  0 o' 

’n  " 

X 

0 

0 4 2 

^2 

= 

0-2  3 

*^2 

+ 

-n(4-8) 

0 2 1, 

■‘^3. 

0 3 -4_ 

.^3. 

-0(0-0) . 

GK*  GC*  W* 


If  the  first  system  shown  in  Equation  1 was  in  progress,  the  element  one  changes 
phase,  then  after  updating,  the  system  should  be  represented  by  Equation  2.  First, 
the  change  in  element  matrices  in  calculated  as  follows: 


'3  4 o' 

'12  0' 

'2  2 0' 

4 3 0 

- 

2 10 

- 

2 2 0 

.0  0 0- 

.0  0 0. 

— 1 

0 

0 

0 

1 

which  is  the  change  for  matrix  ^K*  in  Equation  1,  and 


N> 

00 

0 

__1 

'-4  3 (T 

'6  5 0' 

8 2 0 

- 

3-4  0 

- 

5 6 0 

.0  0 0. 

-0  0 0. 

.0  0 0. 

which  is  the  change  for  matrix  GC*  in  Equation  1. 
Therefore,  Equation  1 becomes 


’3 

2 

0' 

i+1 

'7 

5 

0 

’n ' 

i 

" 0 

2 

4 

2 

as 

5 

-2 

3 

^2 

+ 

-n(2-3) 

_0 

2 

1. 

0 

*1 

-^3 

.-n(o-o). 

(3) 
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Reapplying  the  boundary  condition  t^  “ h*  Equation  3 becomes 


i+1 


O 

O 

1 

1 0 o' 

n" 

0 

0 4 2 

^2 

= 

0-2  3 

^2 

+ 

-n(2-3)-n(2-5) 

.0  2 1, 

0 3 -4. 

.-n(O-O)  -n(O-o) 

0 

4 

2 


■^1 

i+1 

1 

0 

o' 

Vi’ 

0 

^2 

■ 

0 

-2 

3 

^2 

+ 

- n(^~8) 

< 

••^3' 

.0 

3 

■^3. 

n (0-0) . 

GK* 


CC* 


TT 


which  is  exactly  the  relation  stated  in  Equation  2. 

In  conclusion,  if  an  element  changes  phase,  the  following  procedure  is 

used : 


(1)  calculate  the  change  in  thermal  parameters,  i.e.  parameter  (new)  - 
parameter  (old) . 

(2)  using  the  change  in  parameters  as  the  thermal  parameter,  rederive 
the  element  matrices,  modify  the  results  per  the  Crank-Nlcolson 
modifications,  and  add  to  the  global  matrix  system. 

(3)  reapply  the  boundary  conditions. 

(4)  update  element  phase  flagging  whether  frozen  or  thawed. 


The  above  procedure,  however,  does  not  consider  the  changes  imposed  on  GJC* 
and  05*  by  the  insertion  of  boundary  conditions,  which  occurs  prior  to  the  time 
advancement  routine. 

If  the  element  which  changes  phase  is  Involved  with  columns  a and  3,  then  we 
must  also  modify  ^ for  the  new  thermal  parameters. 
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i 

1 

J 


r 


Intuitively,  columns  a and  0 represent  the  columns  of  summed  matrices. 
Hence,  by  the  summing  of  matrices  QC-£)*  and  (X“Z)*  to  G^K*  and  (JC*,  (after 
the  boundary  conditions  were  already  Inserted  in  G^K*  and  QP*) , we  see  that  only 
a column  of  (X“C)*  or  (X“£)*  values  would  show  in  a and  P . If  we  were  to 
reapply  the  boundary  condition  Insertion  scheme,  the  vector  "ij  would  be  properly 
adjusted. 
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TWO-DirvliNGIONAL  TKIANGULAK  iLni  jiNT  "i" 


NTHAV/0(i)*  6 WTHAV.’Dd)  ■ 5 NTriAWi.(i)  ■ 


O thawed  node 
• frozen  node 


(i)  if-!  ad  juste  1 as  follows: 

a)  if  node  A freezes,  add  -1  to  NTHAWD(i) 

b)  if  node  A thaws,  add  ♦ 1 to  NTHAV(’J(i) 

where  "i"  represents  all  elements  associated  to  node  A 

Thus,  if  a node  A changes  phase,  each  element  "i"  containing 
node  A is  determined,  NTHAWD  (i)  is  updated,  and 

a)  if  NTHA'/V7j(  i )s  0 , and  the  element  "i"  is  already 

frozen:  no  phase  change. 

b)  if  NTI{AWr)(i)s  h , and  the  element  "i"  is  already  thawed: 
element  "i"  freeze 

c)  if  NTHAWD(i)“  6 , and  the  element  "i"  is  already  frozen: 

[ element  "i"  thaws. 

1 

d)  if  NTHAWJ(i)e  6 , and  the  element  "i"  is  already  thawed: 

no  phase  change. 


Hence,  if  NTHA\VD(i)  is  between  0 and  6 , there  is  no  need  to 

even  test  for  phase  change.  (Hor  three-dimensional  problems,  the 
critical  numbers  are  HTHAWD(i)  * 0 and  NTHAWL)(i)  ■ 10). 

FIGUKd  C-1 


i 

i 

\ 


111 


APPaNLIX  D 


'"section 

1 

2 

3 

4 

5 


Urj.v.,’Z  i.iAHUAL 

TABL..  OF  GONTjiWTS 
jescription 

3~U  AWIBOTKOPIC  TKANSIKNT  H£AT  CON.^UCTIOF 

WiVjj.'L 

2- J  ANISOTROPIC  TKAN3 ISN'T  H£AT  CON..UCTIC N 
IlOOiiL 

3- .J  ISOTROPIC  H3AT  CONJUCTION  MO.^SL  iVITH 
loOTHFRivlAL  PHASE  CHANGE 

2-D  ISOTROPIC  HEAT  CONDUCTION  MOuEL  WITH 
loOTHiiRMAL  PHASE  CHANGE. 

NUiViBEHING  SEQUENCES 
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— — . — ^ 

o£CTION  1 


FIRST  CONTROL  CARD  (12) 

cc.  2 KOJE:  DIMENSION  OF  PROBLEM  (2) 

SECOND  CONTROL  CARD  (11) 

cc.  1 KODE  1:  "2''  IS  FOR  TRANSIENT  HEAT  CONDUCTION  MOuEL 

(ANISOTROPIC) 


FIRST  DATA  CARD  SET  (6l4,  2F10.5) 


cc . 

1-4 

cc . 

5-8 

cc . 

9-12 

cc. 

13-16 

cc. 

17-20 

cc . 

21-24 

cc . 

25-34 

cc . 

35-44 

NBAND:  BANDWUTH  OF  GLOBAL  MATRICES 

NNODES:  NUMBER  OF  NODES  IN  PROBLEM 

NELE:  TOTAL  NUMBER  OF  ELEi^ENTS 

NBRICK:  NUMBER  OF  BRICK  ELEMENTS 

NTETRA:  NUMBER  OF  TETRAHEDRON  ELEMENTS 

NUMBC:  NUiVBER  OF  NODES  WITH  SPECIFIED  TEMPERATUR.1S 

THETA:  TIME  STEP,  EXPRESSED  IN  HOURS 

DAYS:  .JURATION  OF  TEST,  EXPRESSED  IN  HOURS 


SECOND  DATA  CARD  SET  (6F11.5) 

cc.  1-66  BC(I):  TEMFERAITIRES  OF  BOUNuARY  CONDITIONS,  EX- 

PRESSED IN  °F.  , READ  IN  NUMIERICAL  ORDER-OF- 
MAGNITUDE  OF  CORRESPONDING  GLOBAL  NODE  NUT-IBERS. 


THIRD  DATA  CARD  SET(lll6) 

cc.  1-66  NBC(I):  GLOBAL  NODE  NUIffiERS  OF  BOUNDARY  CONDITION, 

REAJ  IN  NUMBERTCAL  ORDER-CF-FDVGNITUjE. 


113 


FOURTH  JATA  CARj  SjlT  ( 12I6  ./,  12I6  , / , 3I6  ) IBRIC.K  (I.J) 

{Oi.IlT  IF  THFRi:  ARii  NO  FRICR  iiLFMBNTS , I.B,  I'lBiaCK*  0) 

INPUT  BRICK  NO'JAL  :3i£QUi2NCF,  IN  i^LBMENT  OR^E'R  OF  i.lAGNTTUufi. 

NOTE:  Bi;lCK  NO'jAL  INPUT  MUET  BE  IN  "oTANDArE  NOIjAL  UiiQUENCE.  " 

JATA  INPUT  I : PUNCHEiJ  ON  THREE  CArjJ  FOR  c-ACK  BHIC;:  ELEiviENT  INTO 
TIIE  MATRIX  IBRlCK  (I.J)  AS  FOLLOWS: 

car;.'  ONE;  (1216) 

cc.  1-72  FIRST  TWELVi'i  NOEES  IN  SEQUENCE. 

CARD  TWO;  (1216) 

cc.  1-72  NEXT  TWELVE  NONES  IN  SEQUENCE. 

CARD  THREE;  (3l6) 

cc.  1-12  LAST  TWO  NODES  IN  SEQUENCE 

cc.  13-18  ELEMENT  NUI..i-;ER 

(NOTE:  REPEAT  THIS  PROCEDURE  FOR  EACH  BRICK  ELEMENT.) 

FIFTH  DATA  CARD  SST(  6F10 . 4 . AFIO  .4 ) BRICK  (I.J) 

(OMIT  IF  THERE  ARE  NO  BRICK  ELEMENTS , I.E.,  IF  NBHICa*  0) 

INPUT  BRICK  FARAl.-ETER  DATA  FOR  EACH  BRICK  INTO  lATRIX 
BRICK  (I,J)  AS  FOLLOW/ 3: 

CARD  ONE;  (6F10.4) 

cc.  1-in  X(l)  : GLOBAL  X-COORDINATES  OF  BRICK  ELEI.ENT  NODE 

CORRESPONDING  TO  STANDARD  BRICK  NODE  NUIBER  ONr,  (FT) 
cc.  11-20  Y(l):  GLOBAL  Y-COORDINATE  OF  SAID  NODE.  (FT) 

cc.  21-30  Z(l):  GLOBAL  Z-COORDINATE  OF  SAlJ  NODE  (FT) 

cc.  31-40  DELX:  X-JIKENSION  OF  BRICK  (FT) 

cc.  41-50  JELY:  Y-DIMENSION  OF  BRICK  (FT) 

cc.  51-60  DELZ:  Z-DIMENSION  OF  BRICK  (FT) 
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CARD  TWO;  (4F10.4) 

XKX;  THERMAL  COfrjUCTlVITY  IN  X DIRECTION 
(DTU/HR.  FT.  - °F) 

XKY:  THERML  CONDUCTIVITY  IN  DIRECTION  OF  Y 
(BTU/HR  -ft.-  °F) 

XKZ:  THERML  CONDUCTIVITY  IN  Z DIRECTION 
(BTU/HR  -FT.-  °F) 

THER:  PRODUCT  OF  ELEFffiNT'S  DENSITY  (LBW/FT^  ) 

AND  SPECIFIC  HEAT  (BTU/LBM  - °F) 

(NOTE:  REPEAT  THIS  PROCEDURE  FOR  EACH  BRICK  iEHENT. ) 

SIXTH  DATA  CARD  SET  (lll6)  ITETRA  (I,J) 

(OMIT  IF  THERE  ARE  NO  TETRAHEDRON  ELElvlENTS : I.E.  IF  NTETRA*  0) 

INPUT  TETRAHr.DRON  NODAL  SEQUENCE  AND  TETRAHEDRON  ELEI.lENT 
NUI.3ER.  NOTE:  NODAL  INPUT  MUST  BE  IN  "STANDARD  NODAL  SEQUENCc" 

CARD  ONE:  (lll6) 

cc.  l-6o  ten  nodes  in  above  SEQUENCE 

cc.  61-66  ELEMENT  NUI.3ER. 

(NOTE:  REPEAT  THIS  PROCEDURE  FOR  EACH 'ETRAHEDRON  ELEMENT.) 

SEVENTH  DATA  CARD  SET  ( 6F10 . 4 ,/6F10 .4 , /.4F10 .4 , TETRA  (I,J)) 

(OMIT  IF  THERE  ARE  NO  TETRAHEDRON  ELEF.JiNTS:  I.E.  IF  NTETRA*  0) 

INPUT  TETRAHEDRON  PARAMETER  DATA  FOR  EACH  TETRAHEDRON  ELEMENT 
INTO  i-ATRIX  (I.J)  AS  FOLLOWS: 

CARD  ONE:  (6F10.4) 

cc.  l-60  X(l),  Y(l),  Z(l),  X(2).  Y(2).  Z(2) 
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cc.  1-10 

cc.  11-20 

cc.  21-30 

cc.  31-40 


CAH’)  TWO;  (6F10.4) 


CC  . 

l-6o 

X(.3) 

CAR  J 

THREE 

^4 FI 0.4) 

cc . 

c 

1 

XKX: 

cc . 

11-pc, 

XXY; 

cc . 

21-30 

XXZ: 

cc . 

31-40 

THER 

HEAT 

(NOTE 

: REPEAT  ^ 

X(.3).  Y(3).  2(3).  X(4),  Y(4),  Z(4), 


nOr_,> 


Y 

Z 


ZIGHTH  DATA  CARD  SiCT  (6F11.5)  BxiGIN  (I) 

cc.  1-66  FIRST  SIX  I/ilTIAL  TEMPZRATURZS.  (°F)  . INPUT  MUST 

BE  IN  NODAL  ORDER  OF  lAGNITUDE.  INPUT  ALL  INITIAL 
TEMPERATURE (INCLUDING  BOUNDARY  CONuITIONS) . 

FOR  !■  1 TO  NNOJES. 
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..C.CTION  2 


FIRST  CONTROL  CARD  ( 12) 

cc,  2 KOJii:  Uir.iiNSION  OF  PHOBLii;'...  (2) 

SECOND  CONTROL  CARD  (II) 

cc.  1 KODiil:  ''2"  IS  FOR  TKANSIaNT  HFAT  CONjUCTION  i-iOjFL 

(ANISOTROPIC) 

FIRST  r;ATA  CARu  SET  (6F11.5) 

cc.  1-4  NBANO:  BANDWIDTH  OF  GLOBAL  l-ATRICfiS 

cc.  5-8  NNODFS:  NUf.fflER  OF  NODES  IN  PROBLEi'.* 

cc.  0-12  NELS;  TOTAL  NUl-iBSR  OF  ELE^iiNTS 

cc.  13-16  NBRICK;  NUI-BER  OF  BRICK  ELEkENTS 

cc.  17-20  NTETRA:  NUl^BER  OF  TETRAHEDRON  EL£fr:ENTS 

cc.  21-24  NUr.BC:  NUIBER  OF  NODES  WITH  SPECIFIED  TEMPERATURES 

cc.  25-34  THETA:  TIk-E  STEP,  EXPRESSED  IN  HOURS 

cc.  35-44  SAYS:  DURATION  OF  TEST,  EXPRESSED  IN  HOURS 

second  data  card  set  (6F11.5) 

cc.  1-66  BC(I):  TEMPERATURES  OF  BOUNDARY  CONDITIONS,  EX- 

PRESSED IN  °F. , READ  IN  NUkERICAL  ORDER  - OF  - 
kiAGNITUDE  OF  CORRESPONDING  GLOBAL  NODE  NUI.3ERS. 

THIRD  DATA  CARD  SET  (III6) 

cc.  1-66  NBC(I):  GLOBAL  NODE  NUMBERS  OF  BOUNDARY  CONDITION, 

READ  IN  NUMERICAL  ORDER-OF-IAGNITUDE. 
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F'OURTH  jATA  CAKj  SET  (10l6)  IBRICK  (I.J) 


(Oi-IT  IP  THPRR  AHK  NO  B.vICK  BLB'i.jTNTS , I.ii.  iraRlCX-  0) 

INPUT  BRICK  NOJAL  GiiQUKNCK,  IN  h'LRl.JiNT  ORJEH  OF  IAGNITUji:,. 

NOTE:  BRICK  NOOAL  INPUT  iVUGT  BE  IN  "STANJARj  NCjAL  SEQUENCE". 

cc.  1-54  NINE  NOEES  OF  Si^-QUENCE 

cc.  55-60  ELEKIENT  NUKBER 

FIFTH  JATA  CARJ  SET  (7F10.4)  BRICK  (I,J) 

(OMT  IF  THERE  ARE  NO  BRICK  ELEl-liNTS,  I.E.  NBRICK*  O) 

INPUT  BRICK  PARAMETER  DATA  FOR  EACH  BRICK, 
cc,  I-IO  X(l):  GLOBAL  X-COORDINATE  OF  BRICK  ELEMENT  NO  IK 

CORRESPONDING  TO  STANDARD  BRICK  NODE  NUMBER  ONE  (FT), 
cc.  11-20  Y(l):  GLOBAL  Y-COOHDINATE  OF  SAID  NODE  (FT), 

cc.  21-30  DELX:  X-DI1V.ENSI0N  OF  BRICK  (FT), 

cc.  31-40  JELY;  Y-DIMENSICN  OF  BRICK  (FT). 

cc.  41-50  XKX:  THERi/iAL  CONDUCTIVITY  IN  X-dIRECTION  (BTU/HR  . FT-°F  ) 

cc.  .51-60  XKY:  THERML  CONDUCTIVITY  IN  Y- DIRECTION ( BTU/HR . FT-°F ) 

cc,  6l-?0  THER:  PRODUCT  OF  ELEMENT'S  DENSITY(LBM/FT^ ) ANd 
SPECIFIC  HEAT  (ETU/L3i..-°F ) 

(NOTE:  REPEAT  THIS  PROCEDURE  FOR  EACH  BRICK  ELEI-o^iNT,  ) 

■SIXTH  DATA  CARJ  SET  (7l6)  ITETRA  (I.J) 

(OMT  IF  THERE  ARE  NO  TRIANGL  fi  ELEHuENTE  ; I .E.  IF  NTETRA  • 0 ) 

INPUT  TRIANGLE  NODAL  SEQUENCE  AND  TRIANGLE  ELEI.jiNT  N'UlBER. 

NOTE:  NOJAL  INPUT  MUST  BE  IN  "STANDARD  NODAL  SEQUENCE." 
cc.  1-36  SIX  NODES  OF  SEQUENCE 

cc,  37-42  TRIANGLE  ELEMENT  NUivBER 

(NOTE:  REPEAT  THIS  PROCEDURE  FOR  iACH  ELEi.iENT.  ) 
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5i£VSNTH  IJATA  CAR'J  SET  (6F10.4,  3F10.4)  TiiTrtA  (I,J) 

(Oi-.IT  IF  THEHii  ARE  NO  TRIANGLE  ELEKENTS , I.E.  IF  NTETRA«  0 ) 
INPUT  TRIANGLE  PARAi.ETER  .jATA  FOR  EACH  TRIANGLE  ELE1..ENT. 
CAR.)  ONE  (6F10.4) 


CC  . 

l-60 

X(l),  Y(l),  X(2),  Y(2).  X(3),  Y(3) 

CARD 

T'.VO 

(3F10.4) 

cc . 

1-10 

XKX:  THERi/Ai.  CONDUCTIVITY  IN  X-dIRECTION 

(BTU/HR-FT-°P) 

cc . 

11-20 

XKY:  THERJ.Ae  COifcUCTIVITY  IN  Y-DIRECTION 

(BTU/HR-FT-°P) 

cc . 

21-30 

THER:  PRODUCT  OF  DENSITY  (LBk/FT^)  AND  SPECIFIC 
HEAT  (BTU/L31.1-°F) 

(NOTE:  REPEAT  FOR  EACH  TRIANGLE  ELEP^'NT) 

EIGHTH  JATA  CARD  SET  (6F11.5)  BEGIN  (I) 

cc.  1-66  FIRST  SIX  INITIAL  TEiviPERATURES . (°F).  INPUT  ivUST 

3E  IN  KOOAL  OR_'ER  OF  PAGNITUJE.  INPUT  ALL  INITIAL 
TEMPERATURES,  (INCLUJIWG  BOUNDARY  CONDITIONS), 

FOR  I = 1 TO  NNOuES. 
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FOUTtTH  UATA  CARj  YE?  ( 1216  ,/,  1216  ,/3l6 ) IBRICK  (l.J) 

(OfaT  IF  THaRj:  ARR  no  brick  iliifvkJNTJ,  I.ii.  NBRICK  • 0) 

INPUT  BRICK  NOPAL  JBQUBNCi,  IN  BLBi.lBNT  OR.>iiH  OF  ..lAGNITUjxi. 

NOTE:  BRICK  NO  .AL  INPUT  i-lU  T Bii  IN  "oTAN-^ARL  NODAL  .SEQUENCE". 

DATA  INPUT  I : PUNCHE  J ON  THREE  CAR  j3  FOR  .-iACH  BRICK  ELEMENT  INTO 
THE  MATRIX  IBRICK  (1,J)  A:.  FOLLOWE: 

CARj  ONE:  (12l6) 

cc.  l-?2  FIRST  TrtE,LV';  NOjES  IN  oiiQUENCE. 

CARJ  TWO:  (12l6) 

cc.  1-72  NEXT  TV.'jiLVE  NODES  IN  SEQUENCE. 

CAR^  THREE:  (316^ 

cc.  1-12  LAST  TWO  NO  mES  IN  SEQUENCE 

cc.  13-18  ELEilENT  NUMBER 

(NOTE:  REPEAT  THIS  PROCEDURE  FOR  EACH  BRICK  ELEMi:.NT.  ) 

FIFTH  DATA  CARj  SET  (7F10.4,  6F10.4)  BRICK  (I.J) 

(0;.;IT  IF  THERE  ARE  NO  BRICK  ELEIte'NTS,  I.E.  NBRICK*  0) 

INPUT  BRICK  PARAMETER  DATA  FOR  EACH  BRICK  AS  FOLLOWS; 

CARD  ONE;  (7F10.4) 

cc.  1-10  X(l):  GLOBAL  X-COORJINATE  OF  BRICK  ELEi-iENT  NODE 

'corresponding  to  STANjARD  BRICK  NOdE  NU.>BeR  ONE(FT). 
cc.  11-20  Y(l):  GLOBAL  Y-COORJINATE  OF  SAID  NODE  (FT;, 

cc.  21-30  Z(l):  GLOBAL  Z-COORDINATE  OF  SAID  NOjE  (FT), 

cc.  31-40  DELX:  X-DIRENSION  OF  BRICK  (FT), 

cc.  41-50  )ELY:  Y-Jll‘ENSION  OF  BRICK  (FT.) 

cc.  51-60  DELZ:  Z-Jliv'.ENSION  OF  BRICK  (FT), 

cc.  61-70  FROZEN  THErtMiAL  CONDUCTIVITY(BTU/HH . FT-°F ) 
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^IRST  CONTROL  CARL^  (12) 

cc.  2 i<OS£;  Di:.:£NoION  OF  PROBL£i.*  (2) 

SECOrrj  CONTROL  CAR  J (II) 

cc.  1 iCODE  l!  "0"  13  FOR  ISOTROPIC  HEAT  CONJUCTIOH  MOjEL 

WITH  I SOTHERiVAL  PHASE  CHANGE. 


FUST  JATA  CARD 


cc . 

1-4 

cc . 

5-8 

cc . 

9-12 

cc . 

13-16 

cc . 

17-20 

cc . 

21-24 

cc . 

25-34 

cc . 

35-44 

SET  (6l4,  2F10.5) 

NBANl):  BANDWIDTH  OF  GLOBAL  IVIATRICES 

NNOjES:  NULiBEH  OF  NODES  IN  PROBLEM 

NELS:  TOTAL  NUlvBER  OF  ELEl.  ENTS 

N3RICK;  NUMBER  OF  BRICK  ELEMENTS 

NTETRA:  NUMBER  OF  TETRAHEDRON  ELEI.JiNTS 

NU[®C:  NUMER  OF  NODES  WITH  SPECIFIED  TEMPERATURES 

THETA:  TIME  STEP,  EXPRESSED  IN  HOURS 

DAYS:  DURATION  OF  TEST,  EXPRESSED  IN  HOURS 


SECOND  DATA  CARD  SET  (6F11.5) 

cc.  1-66  BC(I):  TEkPERATURES  OF  BOUNDARY  CONDITIONS,  EX- 

PRESSED IN  °F.,  READ  IN  NUlvlERICAL  ORDER-OF- 
MAGNITUDE  OF  CORRESPONDING  GLOBAL  NODE  NW^ERS. 


THIRD  DATA  CARD  SET  (lll6) 

cc.  1-66  NBC(I):  GLOBAL  NODE  NUlfflERS  OF  BOUNDARY  CONDITION, 

READ  IN  NUMERICAL  0RDER-0F-I«1AGNITUDE. 
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CA'A  ) TWO  (6P10.4) 


CC  . 

1-10 

UNFROZEN  THr.Ri’AL  CONDUCTIVITY 

cc . 

11-20 

FROZEN  HEAT  CAPACITY 

cc . 

21-30 

UNF.IOZEN  HEAT  CAPACITY 

cc . 

31-40 

FROZEN  ELc-MENT  DENSITY 

cc . 

41-50 

UNFROZEN  element  jENSITY. 

cc . 

51-60 

(NOTE: 

LATENT  HEAT  PER  UNIT  WEIGHT. 

REPEAT  THIS  PROCEDURE  FOR  EACH  BRICK  ELEf-lENT 

-JIXTH  JATA  CkAO  idT  (111'"')  ITiiTRA  (I,J) 

(OMIT  IF  THliRF  ARi)  NO  TiiTRAHFJRON  ; I.ii.,  IF  NTii'TRA"  0 ) 

INPUT  TKTRAHFiJKON  NODAL  o^QUiiNCii  AND  TETRAHEDRON  ELEiVlENT 
NUi.'lBER.  NOTE;  NOjAL  INPUT  MUST  BE  IN  "STANDARD  NODAL  SEQUENCE" 
CARD  ONE;  (lll6) 

cc.  l-6n  TEN  NODES  IiJ  ABOVE  SEQUENCE 

cc.  61-66  ELEMENT  NU;3ER. 

(NOTE:  REPEAT  THIS  PROCEDURE  FOR  EACH  TETRAHEDRON  ELEMENT.) 


SEVENTH  JATA  CARj  SET  (6F10.4,  6F10.4,  6F10.4)  TETRA  (I,J) 

(or  IT  IF  THERE  ARE  NO  TETiDVHEDRON  ELEMENTS,  I.E.,  IF  NTETKA«0) 
INPUT  TETRAHEDRON  PARAMETER  DATA  FOR  EACH  TETRAHEDRON  ELEI.ENT. 


CARD 

ONE  (6F10.4) 

cc . 

1-60 

X(l),  Y(l),  Z(l).  X(2), 

Y(2),  Z(2) 

CARD 

TWO 

(6F10.4) 

cc . 

1-60 

X(3),  Y(3).  2(3).  X(4), 

Y(4).  Z(4) 

CARD 

THREE 

(7F10,4) 

cc . 

1-10 

FROZEN  ELEMiiNT  CONDUCTIVITY  (BTU/HR-FT-°F) 

cc . 

11-20 

UNFROZEN  ELEMENT  CONDUCTIVITY 

cc. 

21-30 

FROZEN  HEAT  CAPACITY 

J 
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UNFROZEN  HEAT  CAPACITY 


cc.  31-40 
cc.  41-50  FROZEN  ELEMENT  JEN3ITY 

cc.  51-60  UNFROZEN  ELivJviENT  ijENEITY 

cc.  61-70  LATENT  HEAT  PER  UNIT  WEIGHT 

(NOTE:  REPEAT  FOR  EACH  TETRAHEDRON  ELEI.IENT) 

EIGHTH  DATA  CAR.)  SET  (6F11.5)  BEGIN  (I) 

cc.  1-66  FIRST  SIX  INITIAL  TEMPERATURES  (°F).  INPUT  MUST 

BE  IN  NODAL  ORwER  OF  MAGNITUDE.  INPUT  ALL  INITIAL 
TEMPERATURES,  (INCLUDING  BOUN.^ARY  CONDITIONS), 

FOR  1*1  TO  NNOdES. 

NINTH  DATA  CARD  SET  (35F2.0) 

INPUT  INITIAL  PHASE  DATA:  "1"  MEANS  ELEMENT  13  THAWED 

"0"  MEANS  ELEiffiNT  13  FROZEN 
cc.  1-70  35  ELEMENT  PHASE  INDICATIONS 

(REPEAT  ABOVE  PROCEDURE  AS  NEEDED) 


f 


..ACTION  4 


FIHr.T  CONTROL  CAR.)  (12) 

cc.  2 KO^S:  Ji;'.aiN.;iCN  OF  PROBLrJM  (2) 

■j^COrP  CONTROL  CAR..)  (II) 

cc.  1 KO.)S  1;  lo  FOR  ANI.30TK0FIC  TKAN.jIFNT  HiiAT 

CONDUCTION  I 0 )iil,  WITH  iDOTHFRi.iAL  PHADii  CHANG4 

FIR.3T  DATA  CARJ  DFT  (6F11.‘^) 

CC.  1-4  NBANj:  BANjWIDTH  OF  GLOBAL  I.ATHIC4.S 

cc.  5-8  NNOdED:  NU..Br;h  OF  NOJ£  J IN  PROBLLivi 

cc.  9-12  NFL4:  TOTAL  NUMBFR  OF  riLLMFNTD 

cc.  13-16  NBRICK:  NU-.IRLH  OF  BRICK  JiLniMciNTD 

cc.  17-20  NTLTRA:  NUkBFH  OF  TETRAHEDRON  ELEMENTS 

cc.  21-24  NUIIBC:  NUIB.-R  OF  NODES  WITH  SPECIFIED  TEMPERATURE.: 

cc.  25-34  THETA:  TIi-.E  .;TEF,  EXPRESSED  IN  HOUR: 

cc.  35-44  jAY::  juration  OF  TEST,  EXPRESSED  IN  HOUR: 

iECOilD  DATA  CARJ  SET  (6F11.5) 

cc.  1-66  BC(I):  TE:.  F..RATURE3  OF  BOUiPARY  CONDITIONS.  EX- 

PRESSED IN  '^F.,  head  in  NUMBERICAL  ORdER-OF- 
ID^GNITUJE  OF  CORRESPONDING  GLOBAL  NODE  NUMBERS. 

THIRD  DATA  CARJ  SET  (III6) 

cc.  1-66  iNffiCd):  GLOBAL  NODE  NUi.lBERo  OF  BOUNDARY  COIPITION, 

READ  IN  NUfP’RICAL  ORDER-OF-MAGNITUdE. 
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FOUliTH  DATA  CA-i.j  Si:;T  (10l6;  IBl'flCK  (I,J) 

(Oi'IT  IF  THEHi  ARE  NO  BRIG,-:  ii'LSMiiNTS,  I.ii'.  1^RICK*0} 

IjNPUT  brick  nodal  BBQUiiNCii,  IN  JiLKMBNT  ORjiK  OF  l.lAGNITU^r; 
NOTE:  BRICK  NOjAL  INPUT  WUeT  BE  IN  "ETANJARD  NOJAL  SEQUENCE", 
cc.  1-54  NINE  NODES  OF  SEQUENCE 

cc.  55-60  ELE^EiNT  NUMBER 


FIFTH  DATA  CARj  SET  (6F10.4,  5F10.4)  BRICK  (I,J) 

(O.MT  IF  THERJi  ARE  NO  BRICK  ELEi/IENTS:  l.E.  IF  NBRICX*0) 

INPUT  BRICK  PARAMETER  DATA  FOR  EACH  BRICK  ELEi-iENT. 

CARj  ONE  (6F10.4) 

cc.  1-10  . X(l);  GLOBAL  X-COORDINATE  OF  BRICK  ELElnENT  NOjE 

CORRESPONJING  TO  STANDARD  BRICK  NODE  NUl.lBER  ONE  (FT) 
cc.  11-20  Y(l);  GLOBAL  Y-COORDINATE  OF  SAID  NODE  (FT), 

cc.  21-30  DELX:  X-DIIv'.EN'SION  OF  BRICK  (FT), 

cc.  31-40  JELY:  Y-UIMENSION  OF  BRICK  (FT) 

CC.  41-50  FROZEN  THERI  AL  CONDUCTIVITY  (BTU/HR-FT-°F) 
cc.  51-60  UNFROZEN  THERI.IAL  CONDUCTIVITY 

CArj  two  (5F10.4) 

cc.  1-10  FROZEN  HEAT  CAPACITY 

cc.  11-20  UNFROZEN  HEAT  CAPACITY 

cc.  21-30  FROZEN  ELEMENT  DENSITY 

cc.  31-40  UNFROZEN  ELEi*;ENT  DENSITY 

cc.  41-50  LATENT  HEAT  PER  UNIT  WEIGHT 

(NOTE;  REPEAT  THIS  PROCEDURE  FOR  EACH  BRICK  ELEMENT) 
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■^IXTH  DATA  CAR  J (7l6)  ITi^TkA  (I.J) 


(0..t]T  IF  THEF.i  ARE  NO  T.KlAHGLr.  ELEnlENTo ; I.E.  IF  NTET.-iA«C) 
INPUT  TRIAMGLr,  NODAL  ..EQUENCr,  AND  TKlANGLi  ELEkLENT  NUi.iBEH. 
NOTE:  nodal  input  i.iU-.T  Bli.  IN  ''BTANDARij  NOdAL  SEQUENCr,." 
cc.  1-36  SIX  NODE';  or  SEQUi:,liCE 

cc.  37-42  TRIANGLE  ELi.l.iENT  NUMBER 

NOTE:  REPEAT  THIS  PROCejUIiE  FOR  EACH  ELEMENT.  ) 

SEVENTH  DATA  CARD  BET  (7F10.4,  6F10.4) 

(OriT  IF  THERE  ARE  NO  THIANGLE  ELEMENTS : I.E.  IF  KTeTRA*^'). 
INPUT  TRIANGLE  PARAMETER  DATA  FOR  EACH  TRIANGLE. 

CARP  ONE  (7F10.4) 

cc.  1-60  X{1),  Yd).  X(2j,  Y(2),  X(3),  Y(3) 

cc.  61-70  FROZEN  THERiAL  CONDUCTIVITY  (BTU/HR-FT-°F) 


CARD  TWO  (6F10.4) 

cc.  1-10  UNFROZEN  THERMAL  CONDUCTIVITY 

cc.  11-20  FROZEN  HEAT  CAPACITY 

cc.  21-30  UNFROZEN  HEAT  CAPACITY 

cc.  31-40  FROZEN  ELEILMIT  DENSITY 

cc.  41-50  UNFROZEN  ELErM,NT  DENSITY 

cc.  51-60  LATENT  HEAT  PER  UNIT  WEIGHT 

(NOTE:  REPEAT  FOR  EACif  TRIANGLE  ELEMiENT) 

EIGHTH  DATA  CARP  SET  (6F11.5)  BEGIN  (I) 

cc.  1-66  FIRST  SIX  Ii'ilTIAL  TEMPERATURES.  (°F).  INPUT  MUST 

BE  IN  NODAL  ORdER  OF  MAGNITUDE.  INPUT  ALL  INITIAL 
TEMPERATURE  ( INCLUDING  BOUNDARY  CONDITIONS), 

FOR  !■  1 TO  NNODES. 
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raNTH  ’yATA  CAR'J  SilT  (35F2.r) 


cc . 


II'iPUT  Il'ilTIAL  PHAS^  TjATA:  "1"  i.oiAi'iJ  uLiirfiiiNT  IG  THAWi;D 

"0"  i.ih’AMG  iiLdViaMT  IG  FROZ^iA 
1-70  35  liLEi'.liiNT  PHAii;  IRDICATIONG 

(HFFSAT  ABOVE  PROCEGUME  AG  NEEJED) 
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S;'.CT]Oii  5 

NUMBi:,k3NG  PATTr.Kl'l  (2-D) 

Triangle  elements  are  to  be  numbered  from  1 to  NTciTRA.  The 
BRICK  elements  are  to  be  numbered  in  increments  of  "2"  starting 
from  (NTETRA+2)  to  ( ( NTETRA^  2 ) ♦ 2 ( NBRICK)  ) . This  is  because 
there  are  two  triangle  elements  in  each  brick  element. 

ELEwtENT  NUlvlBi'-RIhG  PATTERN  (3-D) 

Tetrahedron  elements  are  to  be  numbered  form  1 to  NTETRA.  The 
BRICK  elements  are  to  be  numbered  in  increments  of  "5"  starting 
from  (NTETRA^  5)  to  ( ( NTETRA-**5 ) + 5 (NBRICK)).  This  is  due  to 
each  BRICK  element  being  composed  of  5 tetrahedron  elements. 

■ITANDARj  NODAL  SEQUEhCEE  (2  or  3-J) 

Refer  to  Chptrs,  2 and  3 for  illustrations  of  the  nodal  numbering 
input  schemes. 
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